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A  Commodity  Supercomputer  For 
Turbulence-Control  Simulations 

DURIP  Award  #:  F49620-99-1-0184 

Jonathan  B.  Freund  and  John  Kim 

Mechanical  and  Aerospace  Engineering 
University  of  California,  Los  Angeles 


The  objective  of  this  DURIP  grant  was  to  build  a  supercomputer  out  of  commodity 
components  to  study  turbulence  and  its  control.  This  objective  has  been  successfully 
met.  This  report  summarizes  the  final  configuration  of  the  system  that  we  built  and 
documents  the  results  of  some  simulations  run  on  it,  including  some  ongoing  AFOSR 
sponsored  work.  The  system  continues  to  be  used  at  a  level  near  its  capacity. 


Configuration 

The  final  configuration  of  the  system  is  as  follows 


Nodes 

20  Compaq  XP1000 

Processors 

Alpha  EV6,  500  MHz 

Memory 

lxl  GB, 16x512  B,  3x256  MB 

Disks 

4  GB  SCSI  each  node 

50  GB  SCSI 

Network 

Fast  Ethernet  Switch 

Operating  System 

Linux  (Red  Hat) 

Software 

LAM  MPI  implementation 
Compaq  Fortran 

800  MHz  Pentium 

Mass  Storage  Unit  1 

500  MB  RAM 

472  GB  Disks  (total) 

2x1  GHz  Pentium 

Mass  Storage  Unit  2 

1  GB  RAM 

472  GB  Disks  (total) 

A  comparison  showed  that  after  vendor  discounts  that  the  20  XPlOOO’s  would,  for  the 
same  price,  offer  comparable  performance  to  Intel  Pentiums  and  greater  flexibility  since 
single-processor  jobs  would  run  three  times  as  fast.  For  many  applications  they  are  still 
faster  than  much  newer  1.5GHz  Pentium  processors. 

Applications 

The  system  has  been  used  to  provide  computation  support  to  AFOSR  Grants 
•  Aeroacoustic  Optimization  and  Control,  Freund,  2000-2002 
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•  Real-Time  Feedback  Control  of  Mixing  in  a  Heated  Jet,  Freund  (via  subcontract), 
1997-2000 

•  Development  of  Robust  Boundary  Layer  Controllers,  Kim,  1997-2000 

•  Robust  Controllers  for  Turbulent  Boundary  Layers,  Kim,  2000-2002 

The  following  are  copies  of  several  publications  resulting  these  projects.  Detailed  techni¬ 
cal  summaries  will  be/have  been  included  in  the  final  reports  for  the  individual  projects. 

Summary  of  publications: 

J.  B.  Freund  &  P.  Moin,  “Jet  mixing  enhancement  by  high  amplitude  fluidic  actuation,” 
with  P.  Moin,  AIAA  J.,  38(10),  1863-1870  (2000). 

P.  Koumoutsakos,  J.  B.  Freund  &  D.  Parekh,  “Evolution  strategies  for  automatic  opti¬ 
mization  of  jet  mixing,”  AIAA  J.  39(5),  967-969  (2001). 

M.  Wei  &  J.  B.  Freund,  “Optimal  Control  of  Free  Shear  Flow  Noise”,  abstract  accepted 
for  the  40th  Aerospace  Sciences  Meeting,  Reno,  Jan  14-17,  2002. 

K.  Lee,  L.  Cortelezzi,  J.  Kim  &  J.  L.  Speyer,  “Application  of  reduced-order  controller 
to  turbulent  flows  for  drag  reduction,”  Phys.  Fluids,  13(5),  2001. 
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Jet  mixing  enhancement  by  high  amplitude  fluidic  actuation 

Jonathan  B.  Freund 

Department  of  Mechanical  and  Aerospace  Engineering 
University  of  California,  Los  Angeles 
Los  Angeles ,  California  90095 
Email:  jfreund@seas.ucla.edu 

Parviz  Moin 

Department  of  Mechanical  Engineering 
Stanford  University 
Stanford,  California  94305 

Abstract 

Recent  experiments  have  shown  that  properly  designed  high  amplitude,  low  mass  flux  pulsed 
slot  jets  blowing  normal  to  a  jet’s  shear  layer  near  the  nozzle  can  significantly  alter  the  jet’s 
development  (Parekh  et  al ,  AIAA  Paper  96-0308).  In  contrast  to  commonly  used  low  amplitude 
forcing,  this  strong  excitation  appears  to  overwhelm  the  turbulence,  having  nearly  the  same  effect 
at  high  and  low  Reynolds  numbers.  It  can  therefore  be  studied  in  detail  by  direct  numerical 
simulation.  In  this  study,  direct  numerical  simulations  of  Mach  0.9,  Reynolds  number  3600 
jets  exhausting  into  quiescent  fluid  are  conducted.  Physically  realistic  slot  jet  actuators  are 
included  in  the  simulation  by  adding  localized  body- “force”  terms  to  the  governing  equations. 
Three  cases  are  considered  in  detail:  a  baseline  unforced  case  and  two  cases  that  are  forced  with 
flapping  modes  at  Strouhal  numbers  0.2  and  0.4.  (St  =  0.4  was  found  to  be  the  most  amplified 
frequency  in  the  unforced  case.)  Forcing  at  either  frequency  causes  the  jet  to  expand  rapidly  in 
the  plane  parallel  with  the  actuators  and  to  contract  in  the  plane  perpendicular  to  the  actuators, 
as  observed  experimentally.  It  is  found  that  the  jet  responds  closer  to  the  nozzle  when  forced 
at  St  =  0.4,  but  forcing  at  St  =  0.2  is  more  effective  at  spreading  the  jet  further  downstream. 
Several  different  measures  of  mixing  (scalar  dissipation,  volume  integrals  of  jet  fluid  mixture 
fraction,  and  point  measurements  of  mixture  fraction)  are  considered,  and  it  is  shown  that  by 
most,  but  not  all,  measures  forcing  at  St  =  0.2  is  the  more  effective  of  the  two  at  mixing. 
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Nomencalature 


a  —  speed  of  sound 

D  =  jet  diameter 

V  =  planar  integral  of  |V£| 

e  =  total  energy 

J  =  volume  integral  of 

M  =  planar  integral  of  £n 

r  —  radial  coordinate 

r0  =  jet  nozzle  radius 

St  =  forcing  Strouhal  number  —  fD/Uj 

T  =  temperature 

Ua  ==  peak  actuator  fluid  velocity 

Uj  —  jet  exit  velocity 

vx  =  axial  velocity 

vr  =  radial  velocity 

vq  =  azimuthal  velocity 

x  =  axial  coordinate 

6  =  azimuthal  coordinate 

p  ==  fluid  density 

^  =  jet  fluid  mixture  fraction 

lo  =  vorticity  magnitude 
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Introduction 


There  are  several  technological  applications  where  enhanced  jet  mixing  can  lead  to  improved  effi¬ 
ciency,  reliability,  or  safety.  For  example,  enhanced  jet  mixing  can  reduce  temperatures  on  in-plume 
aerodynamic  surfaces,  such  as  the  blown  flap  on  a  C-17  aircraft,  and  thus  provide  greater  flexibility 
in  the  choice  of  materials  for  their  construction.  Similarly,  the  mixing  efficiency  of  fuel  jets  in  com¬ 
bustors  is  an  important  factor  in  their  overall  performance,  with  size  and  weight  reductions  possible 
if  mixing  is  improved.  In  the  present  work  we  focus  on  free  jets,  with  the  principle  objective  being 
plume  temperature  reduction. 

In  general,  attempts  to  control  jets  can  be  divided  into  two  categories:  active  and  passive. 
Examples  of  passive  control  are  tabs  located  at  the  nozzle  exit,1,2  crown  shaped  nozzles,3  and 
various  other  tailorings  of  the  nozzle  exit.4-6  This  list  is  far  from  exhaustive.  Passive  control  is 
attractive  because  in  many  cases  it  entails  only  simple  design  modifications,  a  change  in  nozzle 
geometry  for  example.  Also  its  simplicity  typically  makes  the  resulting  hardware  less  subject  to 
failure.  However,  active  control,  where  nozzle  conditions  are  continuously  updated,  has  greater 
flexibility  and  therefore  greater  potential  to  modify  the  jet  flow.  In  this  study  we  analyze  a  recently 
proposed  method  for  active  control  of  high  Reynolds  number  jets.7 

In  the  past,  active  excitation  has  been  used  extensively  to  understand  the  dynamics  of  free  shear 
flows,  particularly  the  dynamics  of  the  largest  turbulent  flow  structures.  Studies  up  to  the  mid 
1980’s  are  summarize  by  Ho  &  Huerre8  and  relatively  more  recent  efforts  are  discussed  by  Koch  et 
al.9  and  Ho  et  al.10  Here,  the  excitation  used  was  typically  low  amplitude,  often  serving  only  to 
seed  instability  waves  in  the  flow  in  order  to  phase  correlate  coherent  structures.  Unfortunately, 
it  seems  that  to  modify  the  flow  significantly  at  high  Reynolds  numbers,  low  amplitude  forcing  is 
ineffective  because  the  applied  perturbations  are  overwhelmed  by  the  turbulence.  Thus,  control  by 
low  amplitude  excitation  is  not  practical  in  many  flows  of  engineering  interest. 
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Figure  1:  Geometry  schematic  showing  the  jet  nozzle  (a)  and  the  actuators  (b). 

Recently,  a  scheme  has  been  developed  to  control  high  Reynolds  number  jets,  such  as  the  exhaust 
flow  from  jet  engines,  by  forcing  with  actuation  velocities  greater  than  the  local  turbulence  intensity. 
This  approach  was  tested  by  Parekh  et  al.7  who  designed  slot  jets  that  blew  normal  to  the  jet  shear 
layer  adjacent  to  the  nozzle  (as  in  figure  1).  When  pulsed  180°  out  of  phase  from  one  another  with 
peak  blowing  at  approximately  one  third  the  jet  velocity  (and  approximately  2  percent  of  the  jet 
mass  flux),  they  excited  large-scale  oscillations  in  the  jet  that  reduced  the  potential  core  length 
by  over  a  factor  of  two.  Recent  test  results  have  shown  mixing  enhancement  using  this  technique 
on  a  full-scale  engine.11  Similar  results  have  been  obtained  using  zero  net  mass  flux  “synthetic 
jets”.7-12-13  It  appears  that  these  active  control  approaches  have  been  more  successful  at  increasing 
mixing  at  high  Reynolds  numbers  than  any  attempts  by  passive  control  approaches.  Obviously  the 
direct  numerical  simulations  used  in  this  effort  are  incapable  of  addressing  high  Reynolds  number 
flows,  but  we  shall  see  that  the  behavior  of  the  simulated  forced  jets  is  similar  to  that  observed  at 
high  Reynolds  numbers,  and  new  insights  are  provided. 

Thus  far,  the  forcing  parameters  in  these  experiments  have  been  picked  beforehand,  as  in  an 
open-loop  control  strategy.  “Closed-loop”  control,  where  the  jet  flow  would  be  continuously 
monitored  and  its  state  used  to  update  the  control  parameters,  may  offer  improved  performance. 
To  implement  this  approach  a  practical  measure  of  the  performance  (an  objective  or  “cost”  function) 
is  needed.  An  objective  of  the  present  effort  is  to  study  metrics  for  mixing  and  provide  a  database 
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for  direct  comparison  of  different  metrics.  The  choice  of  a  metric  for  a  particular  application  will 
depend  not  only  upon  its  relevance  for  the  given  mixing  objective  but  also  on  practical  aspects  of 
its  implementation.  Here  we  concentrate  only  upon  the  metrics  themselves. 

Flow  parameters 

The  focus  of  this  paper  is  a  round  jet  at  Mach  0.9.  The  jet  Reynolds  number,  based  upon  centerline 
flow  conditions  at  the  nozzle  exit,  is  3600,  and  the  stagnation  temperature  of  the  jet  is  constant. 
These  parameters  match  those  studied  experimentally  by  Stromberg  et  a/.14  and  a  baseline,  unforced 
case  has  been  validated  against  their  data.15  Direct  comparison  shows  that  mean  Mach  number 
profiles  (and  the  overall  sound  pressure  levels  on  an  arc  at  60  radii  from  the  nozzle)  agree  to  within  5 
percent.  Unfortunately,  the  precise  nozzle  conditions  were  not  measured,  but  spectra  show  that  the 
initial  jet  shear  layers  were  laminar,  as  expected  at  this  Reynolds  number.  To  model  appropriate 
nozzle  conditions  a  rounded  ‘top-hat’  velocity  profile  was  specified  in  the  present  simulations  (see 
the  Appendix).  Small  amplitude  ( v '  <  O.OlUj)  random  velocity  fluctuations  were  added  to  this  in 
order  to  seed  the  turbulence.  The  consequence  of  not  adding  random  perturbations  was  a  prolonged 
region  of  laminar  flow  near  the  nozzle,  but  the  flow  downstream  was  not  particularly  sensitive  to  the 
nature  of  these  disturbances  provided  that  they  contained  a  range  of  frequencies  and  wavenumbers 
Slot  jets  pictured  in  figure  1  were  used  to  excite  the  flow  in  a  manner  similar  to  the  experiments 
of  Parekh  et  al.7  Each  slot  extended  90°  around  the  jet  and  blew  normal  to  the  shear  layers  just 
downstream  of  the  nozzle.  The  techniques  for  including  these  actuators  in  the  simulations  and 
their  exact  specifications  are  outlined  with  the  numerical  method  in  the  Appendix.  The  individual 
slot  jets  blew  180°  out  of  phase  from  one  another  to  excite  a  flapping  mode  in  the  jet  and  their 
velocity  varied  sinusoidally  between  0  and  0.6 Uj.  The  net  mass-flow  fraction  of  the  actuators  was 
Mact / A/jet  ~  0.035.  Two  forced  jets  where  computed.  The  first  was  forced  at  Strouhal  number 
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Figure  2:  Vorticity  magnitude  contours  for  the  unforced  case  in  an  x-r  plane.  There  axe  40  evenly  spaced  contours 
with  maximum  ur0/Uj  =  8.43. 


St  =  0.2,  which  had  been  found  experimentally  to  be  very  effective  at  spreading  the  jet.7  The  other 
case  was  forced  at  St  =  0.4,  which  was  the  most  amplified  frequency  of  the  unforced  jet,  as  found 
both  by  the  simulations  and  the  cited  experiment.14 


Results 


Visualizations 

Figure  2  shows  a  vorticity  magnitude  visualization  of  the  unforced  case.  The  initial  jet  shear 
layers  are  seen  to  be  laminar.  By  x  =  5r0  instability  waves  appear  which  develop  into  Kelvin- 
Helmholtz  rollers  by  x  =  10ro.  Their  passing  frequency  is  St  =  0.4,  in  accord  with  the  experiments 
of  Stromberg  et  a/.14  who  found  that  St  =  0.44  was  the  peak  Strouhal  number  in  the  early 
development  of  their  jet,  subject  to  ±10%  day-to-day  variation  in  their  facility.  At  the  instant 
shown  the  potential  core  extends  to  approximately  x  =  17r0.  Near  the  end  of  the  potential  core 
a  transition  to  turbulence  occurs  which  is  corroborated  by  Reynolds  stress  statistics  that  will  be 
reported  elsewhere.15 

Figure  3  shows  a  close-up  of  an  actuator  near  its  peak  blowing  condition.  (The  St  =  0.2  case 
is  shown,  but  the  St  =  0.4  case  is  similar.)  There  were  only  8  mesh  points  across  the  modeled 
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Figure  3:  Actuator  flow  impinging  on  the  jet  shear  layer.  Every  other  mesh  point  in  x  and  every  sixth  mesh  point  in 
r  is  shown.  The  solid  rectangle  indicated  the  extent  of  the  actuator.  Light  contours  show  (p  —  Poo)/ PjUj  at  evenly 
spaced  intervals  (min  =  0.0333,  max  =  0.133). 


actuator  in  the  streamwise  direction  and  those  near  the  center  of  the  actuator,  where  the  Gaussian 
axial  velocity  distribution  peaks  (see  the  Appendix),  carry  most  of  the  momentum  flux.  A  region 
of  increased  pressure  just  below  the  actuator  is  also  shown  with  contours  in  the  figure.  The  initial 
effect  of  the  forcing  on  the  jet  is  seen  just  downstream  of  the  actuator  where  the  primary  jet  flow 
is  slightly  deflected.  The  bulk  of  the  fluid  exiting  the  actuator  appears  to  be  turned  downstream 
as  it  encounters  the  jet,  however  a  portion  of  it  is  also  turned  upstream,  giving  the  appearance  of 
a  stagnation  point  flow.  The  pressure  rise  in  this  region  is  also  reminiscent  of  a  stagnation  point 
flow. 

The  result  of  the  forcing  downstream  is  visualized  in  figures  4  and  5  for  the  St  =  0.2  and  0.4 
cases  respectively.  It  is  clear  for  the  St  =  0.2  case  that  the  actuators  excite  the  jet  significantly 
and  that  it  spreads  rapidly  in  the  9  =  -k/2  plane.  In  the  9  =  0  plane,  the  spreading  is  suppressed. 
A  single  mode  appears  to  dominate  the  early  development  of  the  jet  as  evidenced  by  the  large 
structures  visible  before  x  —  10ro.  The  instantaneous  visualization  in  figure  4  shows  the  potential 
core  closing  at  approximately  x  =  9r0.  Downstream  of  this  region,  the  vorticity  magnitude  field 
shows  an  eruption  of  small-scale  turbulence.  A  surprising  feature  of  the  visualization  is  the  apparent 
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Figure  4.  Vorticity  magnitude  contours  for  the  cased  forced  at  St  =  0.2.  Twenty  evenly  spaced  contours  between 
wro/Uj  =  0.6  and  12.0  are  shown. 


symmetry  near  x  —  7 r0  in  the  0  —  0  plane.  This  is  the  characteristic  feature  of  the  large-scale 
coherent  structures  seen  in  the  0  =  7r/2  visualization  as  they  intersect  the  0  =  0  plane. 

When  the  jet  is  forced  at  St  =  0.4  (figure  5),  the  large-scale  structures  are  smaller  but  appear 
earlier,  which  is  not  surprising  because  this  is  the  most  amplified  frequency  in  the  unforced  jet. 
However,  the  downstream  effect  of  the  forcing  is  now  quite  different.  The  structures  disappear 
or  are  obscured  by  small  scales  almost  immediately  and  the  jet  spreads  less  in  the  plane  of  the 
actuators  (0  =  7t/2).  Though  the  shear  layers  appear  thicker  early  in  the  development,  they  slow 
their  spreading  and  merge  only  at  around  x  =  13r0,  beyond  where  they  merge  in  the  St  =  0.2  case. 
The  visualization  in  the  0  =  0  plane  is  similar  to  the  St  =  0.2  case,  with  similar  symmetries  at 
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Figure  5:  Vorticity  magnitude  contours  for  the  case  forced  at  St  =  0.4.  Twenty  evenly  spaced  contours  between 
uir0/Uj  =  0.6  and  12.0  axe  shown. 

small  x  due  to  the  coherence  of  the  excited  structures,  but  a  longer  potential  core. 

Figure  6  shows  a  series  of  instantaneous  visualizations  of  the  jet  fluid  mixture  fraction  as  it 
adjusts  to  forcing  at  St  =  0.2.  The  time  series  starts  from  the  unforced  jet  and  the  actuators 
turn  on  at  t  —  0.  By  t  =  10ro/Uj,  there  is  clear  evidence  of  large  coherent  structures  distorting 
the  scalar  field.  As  expected,  these  travel  at  approximately  0.65f7j,  appearing  at  x  ss  10ro,  the 
horizontal  midpoint  of  the  region  shown,  at  t  ~  1  or()/Uj .  The  scalar  field  takes  considerably  longer 
to  develop  further  downstream  where  decreasing  velocities  slow  advection.  For  computing  statistics, 
the  forced  jets  were  assumed  to  be  fully  developed  only  after  t  =  80ro/Uj.  The  simulations  were 
run  to  approximately  t  =  150ro/Uj. 
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t  =  0ro/Uj  t  =  5r0/Uj 


Figure  6:  Development  of  the  scalar  field  in  the  6  =  7r/2  plane  once  St  =  0.2  forcing  is  turned  on  at  t  =  0.  Black  is 
pure  jet  fluid,  white  is  pure  ambient  fluid. 
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t  =  0ro/Uj  t  =  5r0/Uj 


t  =  10  r0/Uj 


t  =  15  r„/Uj 


t  =  20  rc/Uj 


t  =  25  r0/Uj 


t  =  30  r0/Uj  t  =  35r0/Uj 


Figure  7:  Development  of  the  scalar  field  in  the  6  =  n/2  plane  once  St  =  0.4  forcing  is  turned  on  at  t  =  0.  Black  is 
pure  jet  fluid,  white  is  pure  ambient  fluid. 

Figure  7  shows  the  corresponding  set  of  images  for  the  jet  forced  at  St  =  0.4.  Again  we  see  that 
at  this  forcing  frequency,  large  structures  appear  closer  to  the  nozzle  than  in  the  St  =  0.2  case. 
The  mixed  regions  in  the  jet  shear  layers  thicken  rapidly,  but  a  tongue  of  pure  fluid  persists  along 
the  domain  centerline.  The  flapping  of  the  jet  column  is  also  less  pronounced. 
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Figure  8:  Mean  axial  velocity  ( vx/Uj )  for  the  unforced  jet. 

Mean  flow 


Mean  axial  velocity  ( vx )  provides  a  more  quantitative  measure  of  the  effect  of  the  forcing  on  the 
jet  development.  Approximately  700  fields  spaced  in  time  by  A Tave  =  0.2 r0/Uj  were  averaged  to 
compute  the  mean.  Symmetries  were  exploited  to  increase  the  statistical  sample.  We  see  that  the 
jet  spreads  as  expected  (see  figure  8),  first  slowly  where  the  shear  layers  are  laminar,  and  then  more 
quickly  near  the  end  of  the  potential  core  where  the  flow  goes  through  transition.  If  vx/Uj  =  0.9  is 
used  to  designate  the  end  of  the  potential  core,  in  this  case  the  potential  core  closes  at  x  =  17 r0, 
which  is  further  downstream  than  would  be  expected  for  high  Reynolds  number  jets,  because  the 
shear  layers  are  initially  laminar  and  therefore  spread  slowly. 

Forcing  at  St  =  0.2  dramatically  alters  the  mean  flow  (figure  9).  In  the  plane  perpendicular  to 
the  action  of  the  actuators  ( 9  =  0)  the  jet  at  first  spreads  more  rapidly  than  the  unforced  case,  but 
then  this  is  reversed  starting  at  the  end  of  the  potential  core  {x  =  8r0).  Only  near  x  —  20 r0  does 
the  10  percent  velocity  contour  extend  to  the  same  radial  distance  as  at  x  =  8 r0.  The  contours 
have  an  unusual  thumb  shape  at  r  =  r0,  x  =  8r0  that  is  caused  by  the  large-scale  structures  seen 
in  figure  4.  Based  on  visualizations  ( e.g .  figure  4),  these  structures  first  intersect  the  0  =  0  plane 
near  r  =  0,  and  thus  they  bring  lower  velocity  fluid  into  that  region  before  the  region  near  r  =  r0. 
So  at  this  downstream  location  the  velocity  is  higher  near  r  =  r0  which  causes  the  appearance  of 
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Figure  9:  Mean  axial  velocity  ( vx/Uj )  for  St  —  0.2  forcing. 

the  ‘thumb’.  In  the  9  =  ir/2  plane,  the  jet  is  seen  to  spread  rapidly  starting  near  x  =  8 r0  and 
this  continues  until  the  end  of  the  computational  domain  at  x  =  20ro.  Parekh  et  al.7  showed  very 
similar  results  for  forcing  at  this  same  frequency. 

When  the  jet  is  forced  at  St  =  0.4  (figure  10),  the  mean  flow  is  markedly  different.  In  the 
0  =  0  plane,  the  jet  only  spreads  until  around  x  =  5.5 ra  before  spreading  is  reversed.  It  does  this 
significantly  closer  to  the  nozzle  than  in  the  St  =  0.2  case.  There  is  again  an  unusual  shape  to 
the  contours  at  this  point,  but  at  this  forcing  frequency  it  is  less  pronounced.  Though  spreading 
in  the  9  =  0  plane  is  reversed  at  smaller  x  than  it  was  for  St  =  0.2  forcing,  the  potential  core  is 
now  longer,  extending  to  x  =  13.5 ra.  In  the  9  =  7r/2  plane,  the  jet  spreads  rapidly  starting  at 
around  5.5r0,  but  spreading  slows  downstream  and  the  jet  does  not  grow  as  much  radially  as  in  the 
St  =  0.2  case  (figure  9). 
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t  =  lOrJUj 


t.  =  lHr^/TTs 


Figure  7:  Development  of  tlie  scalar  field  in  the  6  =  ir/2  plane  once  St  =  0.4  forcing  is  turned  on  at  t  =  0.  Black  is 
pure  jet  fluid,  white  is  pure  ambient  fluid. 

Figure  7  shows  the  corresponding  set  of  images  for  the  jet  forced  at  St  =  0.4.  Again  we  see  that 
at  this  forcing  frequency,  large  structures  appear  closer  to  the  nozzle  than  in  the  St  =  0.2  case. 
The  mixed  regions  in  the  jet  shear  layers  thicken  rapidly,  but  a  tongue  of  pure  fluid  persists  along 
the  domain  centerline.  The  flapping  of  the  jet  column  is  also  less  pronounced. 
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Figure  8:  Mean  axial  velocity  ( vx/Uj )  for  the  unforced  jet. 

Mean  flow 

Mean  axial  velocity  (vx)  provides  a  more  quantitative  measure  of  the  effect  of  the  forcing  on  the 
jet  development.  Approximately  700  fields  spaced  in  time  by  Arave  =  0.2ro/Uj  were  averaged  to 
compute  the  mean.  Symmetries  were  exploited  to  increase  the  statistical  sample.  We  see  that  the 
jet  spreads  as  expected  (see  figure  8),  first  slowly  where  the  shear  layers  are  laminar,  and  then  more 
quickly  near  the  end  of  the  potential  core  where  the  flow  goes  through  transition.  If  vx/Uj  =  0.9  is 
used  to  designate  the  end  of  the  potential  core,  in  this  case  the  potential  core  closes  at  x  =  17 r0, 
which  is  further  downstream  than  would  be  expected  for  high  Reynolds  number  jets,  because  the 
shear  layers  are  initially  laminar  and  therefore  spread  slowly. 

Forcing  at  St  =  0.2  dramatically  alters  the  mean  flow  (figure  9).  In  the  plane  perpendicular  to 
the  action  of  the  actuators  (0  =  0)  the  jet  at  first  spreads  more  rapidly  than  the  unforced  case,  but 
then  this  is  reversed  starting  at  the  end  of  the  potential  core  ( x  =  8 r0).  Only  near  x  =  20 r0  does 
the  10  percent  velocity  contour  extend  to  the  same  radial  distance  as  at  x  =  8r0.  The  contours 
have  an  unusual  thumb  shape  at  r  =  r0,  x  =  8 r0  that  is  caused  by  the  large-scale  structures  seen 
in  figure  4.  Based  on  visualizations  ( e.g .  figure  4),  these  structures  first  intersect  the  0  =  0  plane 
near  r  =  0,  and  thus  they  bring  lower  velocity  fluid  into  that  region  before  the  region  near  r  =  r„. 
So  at  this  downstream  location  the  velocity  is  higher  near  r  =  ra  which  causes  the  appearance  of 
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Figure  9:  Mean  axial  velocity  ( vx/Uj )  for  St  =  0.2  forcing. 

the  ‘thumb’.  In  the  9  =  7r/2  plane,  the  jet  is  seen  to  spread  rapidly  starting  near  x  =  8 r0  and 
this  continues  until  the  end  of  the  computational  domain  at  x  =  20 rQ.  Parekh  et  al.7  showed  very 
similar  results  for  forcing  at  this  same  frequency. 

When  the  jet  is  forced  at  St  =  0.4  (figure  10),  the  mean  flow  is  markedly  different.  In  the 
9  =  0  plane,  the  jet  only  spreads  until  around  x  =  5.5r0  before  spreading  is  reversed.  It  does  this 
significantly  closer  to  the  nozzle  than  in  the  St  =  0.2  case.  There  is  again  an  unusual  shape  to 
the  contours  at  this  point,  but  at  this  forcing  frequency  it  is  less  pronounced.  Though  spreading 
in  the  9  =  0  plane  is  reversed  at  smaller  x  than  it  was  for  St  =  0.2  forcing,  the  potential  core  is 
now  longer,  extending  to  x  =  13.5r0.  In  the  9  =  7r/2  plane,  the  jet  spreads  rapidly  starting  at 
around  5.5r0,  but  spreading  slows  downstream  and  the  jet  does  not  grow  as  much  radially  as  in  the 
St  =  0.2  case  (figure  9). 
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Figure  10:  Mean  axial  velocity  (vx/Uj)  for  St  =  0.4  forcing. 

Unsteady  response  of  the  jet 

To  optimize  the  forcing,  by  either  an  open-  or  closed-loop  approach,  it  is  necessary  to  develop 
a  measure  of  its  effectiveness.  Naturally,  the  choice  of  a  definition  for  mixing  effectiveness  will 
depend  upon  the  specific  objective  of  the  application.  Here  we  will  consider  several  metrics  that 
are  of  potential  use  in  temperature  abatement  or  combustion  applications. 

Point  measurements  of  £ 

We  first  consider  point  measures  of  jet  fluid  mixture  fraction  (£  =  1  for  pure  jet  fluid  and  £  =  0 
for  pure  ambient  fluid)  on  the  jet  axis.  If  the  jet  were  hot,  the  concentration  of  jet  fluid  (mixture 
fraction)  would  closely  correspond  to  temperature.  Extensive  centerline  measurements  have  been 
made  in  forced  jets  and  have  been  used  to  estimate  the  mixing  enhancement  of  various  forcings.7 
Figure  11  shows  time  histories  of  (£)  on  the  jet  axis  (r  =  0)  and  at  x  =  5, 10, 15  and  20ro.  The 
angle  braces  ()  indicate  an  average  over  a  single  period  of  the  forcing:  Ta  =  9.9r0/Uj  for  St  =  0.2 
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and  Ta  =  4.9 r0/Uj  for  St  —  0.4.  Despite  this  average,  there  is  still  considerable  oscillation  in  the 
measure  due  to  the  chaotic  nature  of  turbulence.  Averaging  for  longer  periods  would,  of  course, 
smooth  the  profiles,  but  for  closed-loop  control  applications  it  is  important  to  be  able  to  quickly 
measure  the  response  of  the  jet  to  changes  in  the  forcing.  Longer  averages  would  slow  the  response 
of  the  metric  to  the  changes  in  the  forcing. 

We  see  in  figure  11  that  before  forcing  is  initiated  at  t  =  0,  the  flow  is  pure  jet  fluid  at  x  =  5 r0 
and  x  =  10ro.  At  x  =  15 r0  there  is  slight  mixing  with  ambient  fluid  and  at  x  =  20,  the  mixture 
fraction  hovers  around  its  long-time  (if  forcing  were  not  initiated)  mean  £  =  0.65.  As  expected 
from  observations  of  the  potential  core  length,  the  mixture  fraction  at  x  =  5 r0  is  unaffected  by  the 
forcing  at  either  Strouhal  number.  At  St  =  0.2  (figure  11  a),  the  mixture  fraction  at  x  =  10ro  is 
the  first  to  respond  to  the  forcing.  It  initially  decreases  to  (£)  =  0.5,  but  rises  again  and  remains 
near  the  mean  £  =  0.75  for  t  >  60ro/Uj.  The  period  averaged  values  at  x  =  15r0  and  20ro  also 
seem  to  overshoot  initially  before  they  settle  to  hover  around  their  apparent  long-time  mean  values 
of  £  =  0.6  and  £  =  0.3  respectively.  The  small  statistical  sample  size  makes  these  values  and 
the  point  where  they  are  reached  somewhat  imprecise.  Forcing  at  St  =  0.4  also  causes  a  greater 
reduction  in  centerline  mixture  fraction  initially  (figure  lib).  The  curve  at  x  =  10ro  initially  dips, 
but  recovers  to  nearly  its  unforced,  pure  jet  fluid  level  by  t  —  70 r0/Uj.  The  curves  at  x  =  15rc 
and  20ro  hover  around  their  mean  values  of  £  =  0.7  and  £  =  0.4.  It  is  unclear  why  this  case  takes 
longer  to  equilibrate  than  the  lower  frequency  forcing. 

The  initial  overreaction  of  the  jet  to  the  forcing  can  be  explained  qualitatively  with  a  simple 
model  where  the  large-scale  structures  are  assumed  to  be  linear  instability  waves.  Given  this  model, 
turbulent  structures  will  grow,  stabilize,  and  decay  as  the  layer  spreads.  If  the  jet  spreads  slowly, 
as  in  the  unforced  case,  there  is  a  long  region  of  growth  before  decay.  So,  given  a  significant  initial 
forcing  amplitude,  the  structures  can  become  quite  intense  by  a  linear  mechanism.  However,  high 
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Figure  11:  Time  histories  of  scalar  on  the  jet  axis: - x  =  5r0; - x  =  10ro; .  x  =  15r0;  and - x  = 

20ro.  (a)  St  =  0.2;  (b)  St  =  0.4.  Each  curve  shows  a  time  average  over  one  forcing  period. 

amplitude  disturbances  will  also  increase  the  spreading  rate  of  the  layer  and  thereby  reduce  the 
distance  over  which  subsequent  disturbances  can  amplify.  Hence,  by  this  model  the  first  forced 
structure  “sees”  a  slowly  spreading  layer  and  thus  can  grow  more  than  subsequent  disturbances 
that  “see”  a  more  rapidly  spreading  layer,  which  explains  the  observed  overshoot. 

This  model  can  also  explain  the  greater  response  of  the  jet  to  St  =  0.2  forcing  than  forcing  at 
its  natural  frequency  St  =  0.4.  Instability  analysis  shows  that  for  thicker  shear  layers,  as  would  be 
present  in  forced  jets,  the  most  amplified  instability  waves  will  have  a  longer  wavelength  and  lower 
frequency.16  Therefore,  it  is  not  surprising  that  St  =  0.2  is  more  successful.  Unfortunately,  linear 
stability  predictions  can  only  loosely  model  the  quantitative  behavior  of  turbulence.  An  accurate 
quantitative  prediction  using  linear  stability  analysis  is  therefore  not  likely  to  be  successful  in 
conjunction  with  the  high  amplitude  forcing  used  in  the  present  study. 
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Volume  integrals  of  £n 


Volume  integrals  of  can  also  provide  a  measure  of  mixing  effectiveness.  For  n  «  4  to  6  this  will 
provide  a  crude  model,  assuming  £  oc  T,  of  the  infrared  signature  of  the  jet.  Time  histories  of 


where  is  the  physical  domain  (x  <  20 r0,r  <  10ro),  are  shown  in  figure  12.  No  time  averaging 
was  necessary  to  smooth  curves  in  this  case.  We  see  that  the  volume  integrals  of  f 1  and  £2  both 
increase  when  the  jet  is  forced.  For  both  cases,  the  £*  curve  shows  that  there  is  60  percent  more  jet 
fluid  in  the  domain.  Despite  this  increase,  £4  and  £6  decrease  thus  indicating  improved  mixing.  The 
£6  curve  decreases  from  its  unforced  value  J0  to  J  =  0.7  J0  for  the  St  =  0.2  case  and  J  =  0.6^ 
for  the  St  =  0.4  case.  It  is  somewhat  surprising  that  the  St  =  0.4  case  shows  better  mixing  by 
this  measure  because  the  opposite  was  predicted  based  upon  centerline  measurements  (figure  11). 
It  may  also  seem  counter  to  the  visualizations  and  mean  flow  data  which  show  greater  spreading 
for  St  =  0.2  forcing.  Because  (1)  depends  on  the  axial  dimension  of  the  computational  box,  it  is 
therefore  important  to  also  estimate  downstream  mixing  based  upon  the  available  data.  (Note  that 
the  finite  radial  dimension  does  not  affect  the  result  because  £  0  by  r  =  10ro,  the  radial  box 

size.) 

We  can  make  such  an  extrapolation  by  computing  mixedness  as  a  function  of  downstream 
distance.  Figure  13  shows  the  planar  contributions  to  J  as  a  function  of  x  for  t  >  80ro/Uj : 

r2n  rl0ro 

Mx{x)  =  /  /  t^rdrdQ.  (2) 

Jo  Jo 

It  now  becomes  clear  that  the  apparent  advantage  of  St  =  0.4  forcing  is  primarily  a  result  of  the 
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Figure  12:  Integrals  of  with  n  =  1 - ;  n  =  2 - ;  n  =  4 .  ;  n  =  6 - .  (a)  St  =  0.2  forcing;  (b) 

St  =  0.4  forcing. 

jet’s  more  rapid  response  to  the  forcing.  At  the  outflow  boundary  we  see  that  the  St  =  0.2  forced 
jet  is  actually  mixed  better  (by  this  measure)  and  Mx  has  a  steeper  slope  (rate  of  mixing)  than 
the  jet  forced  at  St  =  0.4.  Though  it  is  not  possible  to  make  firm  judgments  about  the  subsequent 
downstream  mixing,  based  on  the  level  and  slope  of  Mx  at  x  =  20ro  it  appears  that  the  St  =  0.2 
case  might  be  better  if  more  downstream  fluid  could  be  included  in  (1).  Both  forced  cases  are 
clearly  better  than  the  unforced  case  also  shown  in  the  figure. 


Figure  13:  Planar  integrals  of  £6  from  (2):  - unforced; - St  =  0.2; .  St  =  0.4. 
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Figure  14:  Planar  integrals  of  £6  from  (3):  - unforced; - St  =  0.2; .  St  —  0.4. 

The  asymmetry  of  the  jets  is  seen  in  figure  14  which  shows  planar  integrals  of  £6  at  constant  9 , 

/20ro  rl0ro 

/  fdrdx.  (3) 

-20 To  Jo 


For  both  forced  cases,  this  metric  peaks  in  the  6  =  ir/2  plane,  the  plane  of  the  actuators.  For 
both  Strouhal  numbers,  this  peak  is  roughly  1.75  times  as  high  as  the  minimum  values.  Somewhat 
surprisingly  these  minima  do  not  occur  at  6  =  0  (see  figure  14)  which  is  perpendicular  to  the 
actuators.  Aie  is  20  percent  higher  at  0  =  0  than  at  its  minima.  Nearly  everywhere  M.g  is 
suppressed  below  its  value  in  the  unforced  case,  more  so  in  the  St  =  0.4  case.  It  should  be  noted 
that  the  details  of  these  results  also  depend  upon  the  box  size.  Because  the  forced  jets  axe  still 
highly  asymmetric  at  the  outflow  boundary,  the  numerical  differences  between  the  peaks  in  valley 
would  increase  if  more  downstream  data  were  available. 
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Figure  15:  Streamwise  mass  flux  from  (4):  - unforced; - St  =  0.2; .  St  -  0.4. 

Streamwise  mass  flux 

The  net  entrainment  of  the  jet  can  be  studied  by  computing  the  streamwise  mass  flux.  Because  vx 
is  negligible  at  r  =  10 r0,  this  is  equivalent  to 

r  rio  rG 

J  pv^rdrdO,  (4) 

which  is  plotted  in  figure  15.  The  mass  flux  in  the  unforced  case  grows  slowly  at  first,  where  the 
shear  layers  are  laminar  and  disturbances  are  small,  and  then  grows  more  rapidly  starting  at  around 
x  =  14r0  where  the  potential  core  closes.  The  mass  fluxes  in  the  forced  cases  both  grow  rapidly 
from  the  start  and  are  over  twice  at  high  by  x  =  20 ra.  Forcing  at  St  =  0.2  is  only  mildly  more 
effective  at  increasing  the  mass  flux  than  forcing  at  St  =  0.4.  It  is  seen  that  the  fluxes  are  slightly 
different  for  the  three  cases  at  x  =  0.  This  is  because  of  the  fluid  added  by  the  actuators  and 
entrainment  caused  directly  by  their  action. 

Scalar  dissipation 

The  rate  of  scalar  dissipation  can  also  provide  an  important  measure  of  mixing  which  is  particularly 
relevant  in  combustion  applications.  In  figure  16,  we  consider  planar  integrals  of  |V£|  as  a  function 
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Figure  16:  Planar  integrals  of  |V£|  from  (5):  - unforced; - St  —  0.2; .  St  =  0.4. 

of  downstream  distance, 


r2ir  rl0ro 

V(x)=  /  |V£|  rdrdd.  (5) 

Jo  Jo 

V(x)  for  the  St  =  0.4  case  starts  to  rise  closer  to  the  nozzle,  but  the  curve  for  St  =  0.2  forcing 
follows  only  2r0  further  downstream  and  becomes  50  percent  greater  by  the  right  side  of  the  domain. 
The  total  dissipation  in  the  computational  domain  is  clearly  larger  for  St  =  0.2  forcing.  Based  on 
the  data  at  x  =  20 r0,  this  trend  is  likely  to  continue  downstream. 


Summary 

Simulations  of  jets  forced  with  high  amplitude  actuation  reproduced  experimental  observations  of 
similarly  forced  jets.  Visualizations  showed  dramatic  effect  of  this  forcing  on  the  jet.  Forcing  at 
St  =  0.2  excited  the  jet  column  into  a  distinct  flapping  mode.  When  forced  at  St  =  0.4,  the  most 
amplified  frequency  in  the  unforced  jet,  the  large  structures  appeared  closer  to  the  nozzle.  However, 
despite  the  rapid  initial  spreading  of  this  jet,  forcing  at  St  =  0.4  was  not  as  effective  at  spreading 
the  jet  and  mixing  it  (by  most  measures)  with  the  ambient  flow  downstream.  Mean  axial  velocities 
showed  that  the  jet  becomes  highly  non-axisymmetric  in  both  forced  cases  and  that  the  potential 
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core  length  was  reduced  more  by  St  =  0.2  than  St  =  0.4  forcing. 

Mixing  was  quantified  by  several  different  metrics.  Point  measurements  of  scalar  concentration 
on  the  jet  axis  showed  that  forcing  at  St  =  0.2  was  more  effective  at  reducing  centerline  scalar 
concentration.  However,  volume  integrals  of  £n  over  the  computational  domain  (for  n  =  4  or  6) 
were  smaller  when  the  jet  was  forced  with  St  =  0.4.  This  was  primarily  due  to  the  faster  response 
of  the  jet  to  forcing  at  this  frequency.  Consideration  of  planar  integrals  of  £6  and  |  V£|  as  a  function 
of  x  suggested  that  forcing  at  St  =  0.2  would  be  more  effective  for  mixing  further  downstream. 
Forcing  at  both  Strouhal  numbers  increased  streamwise  mass  flux  considerably  over  the  unforced 
case,  with  St  =  0.2  forcing  performing  marginally  better  in  this  regard  than  St  =  0.4  forcing. 
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Appendix:  Simulation  techniques 

The  full  compressible  Navier-Stokes  equations  and  an  advection-diffusion  equation  were  solved  in 
cylindrical  coordinates  without  modeling  assumptions.  Full  details  of  the  basic  numerical  algorithm 
are  given  by  Freund  et  al.17  where  the  same  algorithm  was  used  to  study  turbulence  in  compress¬ 
ible  mixing  layers.  Only  a  summary  of  the  method  is  provided  here.  Sixth-order  compact  finite 
differences18  were  used  to  compute  derivatives  in  the  streamwise  (x)  and  radial  (r)  directions  and 
Fourier  spectral  methods  were  used  in  the  azimuthal  direction  (6).  A  fourth-order  Runge-Kutta 
algorithm  was  used  to  advance  the  solution  in  time.  At  the  r  =  0  coordinate  singularity,  the  equa¬ 
tions  were  solved  in  Cartesian  coordinates.  To  maintain  a  reasonable  numerical  time  step  given  the 
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Figure  17:  Mesh  showing  every  tenth  mesh  point  for  the  forced  jet  simulations.  The  thick  rectangle  demaxks  the 
physical  portion  of  the  computational  domain. 

restriction  imposed  by  the  Runge-Kutta  algorithm,  higher  Fourier  modes  in  6  were  not  used  near 
r  =  0.  They  were  omitted  systematically  so  that  the  effective  azimuthal  resolution  remained  nearly 
constant  with  radial  location.  The  computational  mesh  had  440  x  230  x  160  points  in  the  axial, 
radial,  and  azimuthal  directions  respectively.  Mesh  points  were  compressed  in  the  radial  direction 
around  r  =  r0. 

The  boundary  conditions,  the  jet  nozzle,  and  the  actuators  were  all  accounted  for  by  modifying 
the  equations  in  appropriate  regions  of  the  computational  domain.  The  rest  of  the  computational 
domain  that  is  free  of  these  terms  is  referred  to  as  the  “physical  domain”.  The  artificial  terms 
serve  to  model  an  infinite  domain,  the  nozzle,  or  the  actuators  where  needed.  For  example,  non¬ 
reflecting  outflow  and  lateral  (large  r)  boundary  conditions  were  approximated  by  stretching  the 
mesh  beyond  the  physical  domain  and  filtering  the  solution  on  that  stretched  mesh.  A  typical  mesh 
is  shown  in  figure  17.  Similar  boundary  conditions  have  been  used  in  aeroacoustic  computations 
and  additional  documentation  is  available  elsewhere.15,19 

In  a  zone  of  width  ra  at  the  inflow  side  of  the  computational  domain  (see  figure  17),  terms 
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were  added  to  the  equations  that  drive  the  solution  toward  the  desired  nozzle  conditions.  If  the 
compressible  flow  equations  are  represented  by  N( q)  =  0,  where  q  is  a  vector  of  flow  quantities, 
then  the  modified  equation  has  the  form  N( q)  =  -<7(q-qtarget),  where  qtargct  is  the  desired  nozzle 
conditions.  The  added  term  acts  like  a  penalty  function.  The  value  of  a  in  this  zone  was  2.50oo /r0 
and  the  target  axial  velocity  was  a  typical  hyperbolic  tangent  ‘top-hat’  profile: 


target 


1  —  tanh 

12.5(—  )1 

_ 

V^O  rJW 

(6) 


The  density  and  scalar  profiles  were  similarly  specified.  Random  velocity  fluctuations  where  added 
to  this  mean  profile  using  a  similar  body- “force”  method.  This  randomization  was  low  amplitude 
( v '  <  O.Olf/j)  and  served  only  to  seed  the  turbulence  with  broad-banded  noise. 

The  actuators  where  included  by  similar  body-force  methods.  The  actuators  extended  axially 
from  x\  =  0.13ro  to  =  0.57ro  (8  mesh  points)  and  radially  from  r\  =  l.llr0  to  r 2  =  3.5r0 
(101  mesh  points).  In  this  region,  a  =  2.5a00/r0  and  qtarget  =  [p, pvx , pvr , pvg , e, K]targot  = 
[•75  p  00, 0,  -75^00^0) 0,  ea, 0]  with  ca  calculated  from  p  —  .75poQ,  and  the  ambient  pressure. 

The  instantaneous  actuator  target  velocity  was 


Va  =  Uafx(x)fr(r)fg(0)ft(t). 


(7) 


Its  amplitude  was  set  at  Ua  =  Uj,  but  the  actual  peak  vr  at  the  actuator  only  reached  O.GUj  because 
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Figure  18:  The  x-r  shape  of  the  actuation  velocity. 


Figure  19:  The  9  shape  of  the  actuation  velocity. 

the  body  forcing  was  resisted  by  viscosity.  The  spatial  and  temporal  dependencies  of  va  were 


fx(x)  =  e— Mz-°.5(x1+z2))2  (g) 

Mr)  =  e~a^-^2  (9) 

M9)  =  exp[-C(0)8e_CW  +2  8]  +  exp[-(C(0)  -  4)8e_(<(tf)_4)2+2.8j 

ft{t)  =  0.5  [l  +  sin(27rS'tl7,t/n)sgn(sin0)] ,  (11) 


where  ax  =  95.1  ,ar  =  1.34,  and  ((d)  =  86>/2t r  -  2.  The  x-r  shape  of  va  is  plotted  in  figure  18  and 
the  azimuthal  dependence  is  plotted  in  figure  19.  The  ‘sgn’  term  in  (11)  makes  the  actuators  act 
out  of  phase. 
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1  Introduction 


Evolution  strategies  (ES)  are  introduced  for  the  optimization  of  active  control  parameters  for  en¬ 
hancing  jet  mixing.  It  is  shown  that  the  evolution  algorithms  can  identify,  in  an  automated  fashion, 
not  only  previously  known  effective  actuations  but  also  find  good  but  previously  unidentified  pa¬ 
rameters.  In  this  study  simulations  of  model  jets  are  used  to  demonstrate  the  feasibility  of  the 
methods.  Evolution  strategies  are  robust,  highly  parallel  and  portable  algorithms  that  may  be 
most  useful  in  an  experimental  setting  at  realistic  Reynolds  numbers.  Simulations  of  inviscid  in¬ 
compressible  flows  using  vortex  models  as  well  as  direct  numerical  simulations  (DNS)  of  very  low 
Reynolds  number  compressible  flows  are  used  in  this  study  to  evaluate  different  forcing  parameters. 

Our  objective  is  twofold:  (i)  explore  the  possibility  of  ES  to  find  previously  identified  modes  of 
efficient  operation,  and  (ii)  discover  previously  unknown  effective  actuations.  Practical  engineering 
concerns  will  dictate  the  choice  of  actuator  parameters  and  relevant  cost  functions.  In  Section  2 
we  present  a  description  of  the  ES;  in  Section  3  we  present  results  from  the  application  of  these  to 
the  optimization  of  compressible  jets  and  vortex  models.  Section  4  is  a  discussion  of  results  and  an 
outline  for  future  research. 

2  Evolution  Strategies 

Evolution  strategies  are  continuous  parameter  optimization  techniques  based  on  principles  of  evo¬ 
lution  such  as  reproduction,  mutation  and  selection.  We  define  a  vector  in  the  control  parameter 
space  x  =  (x\,X2 as  an  individual  and  a  set  of  such  individuals  as  a  population.  Evolu¬ 
tion  strategies  use  a  fitness  value ,  prescribed  by  F(x)  =  F(xi,X2,  to  identify  the  best 

individual  from  a  population.  We  take  better  individuals  to  have  larger  F  values. 

2.1  Two-membered  Evolution  Strategies 

The  simplest  evolution  strategy  has  a  population  with  two  competing  individuals,  a  two-membered 
strategy .  Evolution  occurs  by  mutation  and  selection ,  the  two  operations  that  Darwin  considered 
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the  most  important  in  the  evolution  of  species.  Each  individual  is  represented  by  a  pair  of  real 
vectors  u  =  u(x,3),  where  3  is  an  M-dimensional  vector  of  standard  deviations. 

Following  Rechenberg,1  the  optimization  algorithm  is  as  follows: 

1.  Initialization:  A  parent  genotype  consisting  of  M-genes  is  specified  initially  as  x°. 

2.  Mutation:  The  parent  of  generation-n  produces  a  descendant  with  slightly  different  genotype. 

The  operation  of  mutation  is  realized  by  modifying  xp  according  to  x%  =  +  J\t(3p), 

where  Af{3)  denotes  an  M-dimensional  vector  of  normal  random  number  with  zero  mean  and 
standard  deviations  3. 

3.  Selection:  The  fittest  individual  according  to  its  F  value  becomes  the  parent  of  the  next 
generation: 

jn+l  =  fa?.  ifF(zJ)>F(f3; 

p  otherwise. 

The  variance  of  the  population  members  is  adjusted  using  the  1/5  success  rule  proposed  by 
Rechenberg:  if  more  than  one  in  five  offsprings  result  in  an  improved  solution,  then  the  variance 
is  increased.1  For  regular  optimization  problems2  the  method  is  known  to  converge  to  a  global 
minimum,  but  the  rate  of  convergence  cannot  be  anticipated.  So  in  finite  time  there  is,  of  course,  no 
guarantee  that  the  global  optimum  has  been  reached,  a  trait  shared  by  all  optimization  techniques. 
Schwefel3  provides  a  complete  description  of  the  algorithm. 

2.2  Parameter  constraints 

In  problems  of  active  flow  control,  engineering  considerations  impose  constraints  on  the  actuation 
parameters.  Such  constraints  are  formulated  as  inequalities,  such  as  Cj(x)  >  0.  In  this  work 
descendants  that  do  not  satisfy  the  constraints  are  treated  as  unsuccessful  mutations. 
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3  Jet  Flows  Optimization 


3.1  Optimized  Excitation  of  Compressible  Jets 

Direct  numerical  simulations  of  the  developing  region  of  compressible  jets  forced  by  slot-jet  fluidic 
actuators  are  used  to  evaluate  the  fitness  of  individuals.  The  compressible  flow  equations  were 
solved  directly  using  a  combination  of  sixth-order  compact  finite  differences,  spectral  methods, 
and  fourth-order  Runge-Kutta  time  advancement.  Further  details  of  the  numerical  algorithm  and 
techniques  for  including  two  slot-jet  actuators  actuators  that  each  span  90°  of  the  jet  (on  opposite 
sides)  just  down  stream  of  the  nozzle  were  reported  elsewhere.4  Despite  the  low  Reynolds  number 
dictated  by  the  computational  expense  (only  Re  =  500  in  this  study)  the  actuators  were  able  to 
induce  the  gross  effects  observed  in  experiments  at  much  higher  Reynolds  numbers.5  A  path  to 
implementation  at  more  relevant  flow  conditions  is  discussed  in  Section  4. 

For  this  study,  only  three  actuation  parameters  were  varied:  the  amplitude,  frequency,  and 
phase.  The  actuation  is  taken  as  a  sum  of  harmonic  waveforms: 


MAX 


1  +  sin  (u-j +  <p?j sgn[cos 6}  ,^, 


(2) 


where  vr  is  the  radial  velocity  at  the  actuator  exit,  Uj  is  the  jet  exit  velocity,  At  are  the  amplitudes, 
Sti  are  the  Strouhal  numbers,  and  (pi  are  the  phases  of  the  different  modes.  The  sgn[cos(0)]  factor 
causes  each  waveform  to  excite  a  flapping  mode  in  the  jet.  Note  that  the  phases,  (pi,  are  the 
relative  phases  of  the  different  modes — the  two  actuators  always  acted  180°  out  of  phase.  The  Ai 
were  constrained  to  be  non-negative  and  the  Sti  were  restricted  to  be  0  <  St  <  0.8.  The  relative 
phases  were  not  constrained. 

The  computational  mesh  was  112  x  42  x  16  in  the  streamwise,  radial,  and  azimuthal  direction 
respectively  and  the  computational  domain  extended  to  16r0  downstream  and  5r0  in  the  radial 
direction,  where  r0  is  the  jet  radius.  A  stretched-mesh  boundary  zone  was  positioned  in  x  >  8 r0 
and  r  >  3.5r0  to  absorb  out-flowing  fluctuations.4  In  each  iteration  of  the  optimization  the  jet  was 
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simulated  starting  from  an  unforced  case  for  several  periods  of  forcing  after  the  passing  of  initial 
transients.  Because  the  flow  is  laminar  and  becomes  quasi-periodic,  this  was  sufficient  to  provide 
a  measure  of  the  long-time  actuator  effectiveness.  In  total  200  iteration  were  made. 

Three  wave  forms  ( N  =  3)  were  used.  The  initial  parameters  and  the  parameters  that  maxi¬ 
mized  F  are  shown  in  Table  1.  The  fitness  function  was  defined  as 


F  = 


v/rdrdddx, 


(3) 


which  with  was  evaluated  using  a  trapezoidal  rule  quadrature  and  increased  by  a  factor  of  10  for 
the  best  parameters. 


Initial 

Best 

Ai/U 

su 

<t>i 

Ai/U 

su 

0.45 

0.50 

0.00 

0.04 

0.33 

0.54 

0.40 

0.20 

0.70 

0.42 

0.17 

0.31 

0.35 

0.50 

1.00 

0.07 

0.45 

1.57 

Table  1:  Initial  and  best  actuation  parameters. 

It  is  interesting  to  note  that  the  evolution  strategy  reduced  the  amplitude  of  two  of  the  wave 
modes  to  a  very  low  level  giving  effectively  the  same  excitation  that  was  shown  to  be  successful 
in  experiments5  and  in  higher  Reynolds  number  simulations  on  larger  meshes.4  This  also  demon¬ 
strates  the  mesh  independence  of  the  computation  and  its  accuracy,  at  least  for  this  measure  of 
mixing.  Notice  also  that  the  maximum  amplitude  stayed  below  half  the  jet  velocity  indicating  that 
a  sinusoidal  profile  was  preferable  to  one  that  was  clipped  by  the  maximum  amplitude  constraint 
in  (2).  Scalar  mixture  fraction  is  visualized  in  figure  1  for  the  initial  and  best  cases.  The  best  case 
clearly  shows  high-amplitude  flapping. 


3.2  Incompressible  Vortex  Model 

In  these  simulations  it  is  assumed  that  the  effects  of  compressibility  and  viscosity  do  not  affect  the 
flow  dynamics  and  the  circular  jet  is  modeled  by  the  combination  of  discrete  vortex  filaments  and  a 
fixed  semi-infinite  cylindrical  sheet  of  vorticity  with  circulation  per  unit  length  7  that  represents  the 
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nozzle.  Helical  excitation  as  used  in  the  experiments  of  Lee  &  Reynolds6  is  modeled  by  rotating 
the  axis  of  the  vortex  cylinder  with  displacement  Ah/r0  about  the  nominal  jet  centerline.  The 
rotation  frequency  is  denoted  by  fh  and  the  axial  frequency  is  defined  as  fa  =  Staj/2ra.  The 
frequency  fa  is  the  rate  at  which  filaments  are  generated  at  the  origin.  The  velocities  induced  by 
each  filament  and  the  vortex  sheet  are  added  to  determine  the  trajectory  of  each  filament.  Sta  sets 
the  time  between  creation  of  new  rings.  The  circulation  of  each  shed  filament  is  7  =  Star/R.  (3 
fixes  the  ratio  of  the  axial  to  orbital  excitation,  and  <f>  sets  their  relative  phase.  Further  details  of 
this  model  and  its  numerical  implementation  were  reported  by  Parekh  et  al.7  which  also  provided 
a  detailed  validation  of  the  method  through  direct  comparison  with  bifurcating  jet  experiments  and 
refinements  of  the  vortex  ring  representations.6 

The  control  parameters  were  constrained  so  that  0  <  Ah  <  1,  0.1  <  Sta  <  1,  0.2  <  /3  <  5, 
0  <  <f>  <  2ir.  The  fitness  function  was  the  average  angle  of  the  ring  trajectories  defined  as  the  angle 
between  the  jet  centerline  and  the  line  that  connects  the  center  of  the  jet  exit  to  the  centroid  of 
the  vortex  ring  nodes.  This  metric  was  evaluated  after  eleven  periods  of  axial  excitation. 

Initially  we  expected  the  algorithm  to  select  a  jet  that  bifurcates  in  a  single  plane  with  values  of 
Sta  and  Ah  that  maximize  the  spreading  angle  as  seen  in  laboratory  jets.6  Instead,  a  jet  flow  was 
found  that  had  never  been  observed  in  previous  experiments  or  calculations.  It  initially  resembles 
a  bifurcating  jet  (Fig.  2),  but  several  diameters  downstream  the  two  branches  of  the  jet  exhibit 
a  secondary  bifurcation  in  which  the  rings  change  direction  by  about  45°.  This  results  in  a  wide 
spreading  angle  as  seen  in  Fig.  2.  Unfortunately,  the  parameters  that  lead  to  this  double  bifurcation 
have  not  been  tested  experimentally  so  the  code  can  not  be  validated  for  this  case  as  is  was  for  the 
other  cases.7  Regardless,  the  possibility  that  ES  can  identify  novel  control  parameters  has  been 
demonstrated. 

4  Summary  and  Conclusions 

Given  the  variety  of  different  forcing  schemes,  cost  functions,  and  flows  to  be  controlled,  it  is 
impossible  to  anticipate  the  best  optimization  scheme  for  any  particular  case,  but  the  present 
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results  clearly  demonstrate  that  ES  can  be  a  valuable  tool  for  jet  mixing  optimization.  They 
identified,  in  an  automated  fashion,  previously  known  flow  controls  and  found  previously  unknown 
parameters  that  further  enhance  spreading.  Their  strength  is  their  portability,  which  is  considered 
particularly  attractive  for  problems  such  as  jet  noise  control  flow  where  the  mechanisms  appear  too 
poorly  understood  to  provide  much  direct  guidance. 

In  computational  studies,  such  as  the  present,  finite  computational  resources  restrict  the  Reynolds 
numbers  that  can  be  addressed.  The  implementation  of  turbulence  models  and  large-eddy  simu¬ 
lation  calculations  is  a  natural  next  step,  and  is  being  carried  out.8  Besides  computations,  ex¬ 
periments  are  appealing  for  many  applications  as  they  can  provide  rapid  answers  at  realistic  flow 
conditions.  Such  a  study  has  been  initiated  by  the  authors  inspired  by  the  present  results. 

The  future  role  of  computations  in  such  optimizations  is  to  address  issues  where  experiments 
are  limited.  For  example,  consider  the  optimization  of  the  physical  shape  of  realistic  actuators. 
It  would  be  difficult  to  design  hardware  with  the  flexibility  to  provide  a  general  shape  for  the 
actuation.  But  in  a  simulation  it  is  straightforward  to  implement  actuators  of  nearly  any  shape, 
and  they  may  be  constrained  so  that  determined  optimal  geometry  is  realizable  in  hardware.  This 
way  the  final  configuration  can  be  built  and  applied. 

References 

[1]  Rechenberg,  I.,  Evolutionsstrategie:  Optimierung  technischer  System  nach  Prinzipien  der  biol- 
ogischen  Evolution ,  Fromman-Holzboog,  Stuttgart,  1973. 

[2]  Michalewicz,  Z .,  Genetic  Algorithms  +  Data  Structures  +  Evolution  Programs ,  Springer- Verlag, 
Berlin,  1996. 

[3]  Schwefel,  H.-P.,  Evolution  and  optimum  seeking ,  Wiley,  1995. 

[4]  Freund,  J.  B.  and  Moin,  P.,  “Jet  mixing  enhancement  by  high  amplitude  fluidic  actuation,” 
AIAA  J.,  Vol.  38,  No.  10,  2000. 


7 


[5]  Parekh,  D.  E.,  Kibens,  V.,  Glezer,  A.,  Wiltse,  J.  M.,  and  Smith,  D.  M.,  “Innovative  jet  flow 
control:  mixing  enhancement  experiments,”  AIAA  Paper  96-0308,  1996. 

[6]  Lee,  M.  and  Reynolds,  W.  C.,  “Bifurcating  and  blooming  jets,”  Report  No.  TF-22,  Department 
of  Mechanical  Engineering,  Stanford  University,  1995. 

[7]  Parekh,  D.,  Leonard,  A.,  and  Reynolds,  W.  C.,  “Bifurcating  jets  at  high  Reynolds  numbers,” 
Report  Nr.  TF-35,  Department  of  Mechanical  Engineering,  Stanford  University,  1988. 

[8]  Hilgers,  A.,  Center  for  Turbulence  Research  Annual  Research  Brief,  chap.  Parameter  optimiza¬ 
tion  in  jet  flow  control,  CTR,  1999,  pp.  179-194. 


8 


x/r0 
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With  the  advent  of  high-speed  processors  and  advanced  numerical  methods,  it  has  re¬ 
cently  become  possible  to  compute  aerodynamic  noise  from  first  principles,  by  solving  the 
compressible  flow  equations  without  modeling  approximation.1  Two-dimensional  mixing  lay¬ 
ers,2  axisymetric  jets,3  and  supersonic4  and  subsonic1  three-dimensional  turbulent  jets  have 
been  computed  in  this  way.  Those  simulations  have  provided  a  considerably  deeper  picture 
of  the  mechanisms  of  sound  generation,  but  as  of  yet  have  not  yielded  a  clear  method  for 
noise  reduction.  In  this  study  we  incorporate  this  simulation  capability  into  a  new  method 
to  analyze  and  control  flow  noise. 

Our  focus  is  jet  noise  and  we  consider  a  two-dimensional  model  of  the  near  nozzle  region. 
The  action  of  an  unsteady  actuator  is  modeled  by  a  body  force  acting  on  the  shear  layer 
near  the  nozzle.  Our  objective  is  to  reduce  radiated  noise  in  a  region  of  the  sound  field  as 
shown  schematically  in  figure  1  (details  in  the  figure  will  be  explained  later).  Using  adjoint- 
based  optimization,5  the  gradient  information  determined  by  a  solution  of  the  adjoint  flow 
equation  is  used  to  optimize  the  forcing  automatically. 

A  high  order  scheme  with  no  dissipation  and  high  resolution  in  wide  spectrum  is  need  to 
catch  the  sound  wave  accurately.  Following  the  successful  application  in  flow  and  acoustic 


1 


simulations  taken  by  Freund4,1  the  sixth  order  Pade  scheme  discussed  by  Lele6  is  used  in 
spatial  differencing,  and  the  forth  order  Runge-Kutta  algorithm  is  used  for  time  advance¬ 
ment.  Both  buffer  zone7  and  non-reflecting  boundary  condition8  are  used  to  avoid  the 
contamination  from  reflected  wave. 

The  methods  are  benchmarked  in  an  anti-sound  problem  as  illustrated  in  Figure  2, 
though  we  stress  that  antisound  will  not  be  used  for  the  mixing  layer  control.  A  noise 
source  sits  in  a  mean  flow  with  Mach  number  0.5.  Our  objective  is  to  minimize 

J  =  f  f  p^dCldt, 

Jo  J n 

where  is  the  vertical  line  in  Figure  2.  Reducing  the  noise  along  this  line  is  equivalent 
to  reducing  cost  function  J.  To  achieve  this  purpose,  a  mass  source  0(x,  t)  of  compact 
support  in  C  is  optimized  as 

<f>k+1  =  <j>k-  akgk , 

where  g  is  the  gradient  information  from  adjoint  approach  and  a  is  a  parameter  of  descent 
which  governs  how  large  an  update  is  made.  In  real  calculation,  this  simple  gradient  update 
is  inefficent  because  of  the  complexity  of  function  </>,  So  the  conjugate  gradient9  is  used  for 
the  control  update,  with  a  computed  at  each  iteration  by  Brent’s  method.9  Finally,  the 
cost  J  is  reduced  by  more  than  95%  in  only  3  iterations  as  shown  in  figure  3. 

As  the  near  nozzle  jet  model,  two-dimensional  mixing  layer  and  its  radiated  sound  field 
are  simulated  (figure  1).  The  flow  above(inside)  the  nozzle  lip  has  Mach  number  0.8,  and  the 
flow  under  the  nozzle  lip  has  zero  velocity.  Vorticity  thickness  is  used  as  characteristic 
length,  and  the  shear  speed  is  used  as  charactersitic  speed.  Then  we  choose  Reynolds 
number  500  by  definition  Re  =  (ui  -  U2)5u/v.  In  this  case,  temperature  profile  is  chosen 
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to  be  uniform.  The  small  square  area  near  the  middle  of  left  boundary  is  the  area  to  put 
actuating  force  to  control  the  acoustic  field.  Freund  et.al’s  work10  shows  us  the  possibility 
of  changing  the  entire  jet  flow  structure  only  by  actuating  within  near  nozzle  region.  The 
line  below  vortices  is  the  target  area  Q  where  cost  function  should  be  minimized. 

The  same  cost  function  J  is  chosen  along  the  line  to  indicate  the  noise  level.  The  cost 
is  reduced  by  about  46%  after  3  iterations  (figure  4),  while  the  maximum  amplitude  of  the 
energy  put  to  control  is  only  about  10%  of  the  energy  of  the  flow.  Figure  5  compares  the  flow 
and  acoustic  field  for  both  the  original  and  optimized  cases  at  the  same  time.  Significant 
sound  reduction  can  be  observed  in  the  target  area.  Futher  improvements  are  expected  as 
the  method  is  refined  and  run  longer. 
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A  reduced-order  linear  feedback  controller  is  designed  and  applied  to  turbulent  channel  flow  for 
drag  reduction.  From  the  linearized  two-dimensional  Navier-Stokes  equations  a  distributed 
feedback  controller,  which  produces  blowing/suction  at  the  wall  based  on  the  measured  turbulent 
streamwise  wall-shear  stress,  is  derived  using  model  reduction  techniques  and  linear- 
quadratic-Gaussian/loop-transfer-recovery  control  synthesis.  The  quadratic  cost  criterion  used  for 
synthesis  is  composed  of  the  streamwise  wall-shear  stress,  which  includes  the  control  effort  of 
blowing/suction.  This  distributed  two-dimensional  controller  developed  from  a  linear  system  theory 
is  shown  to  reduce  the  skin  friction  by  10%  in  direct  numerical  simulations  of  a  low-Reynolds 
number  turbulent  nonlinear  channel  flow.  Spanwise  shear-stress  variation,  not  captured  by  the 
distributed  two-dimensional  controller,  is  suppressed  by  augmentation  of  a  simple  spanwise  ad  hoc 
control  scheme.  This  augmented  three-dimensional  controller,  which  requires  only  the  turbulent 
streamwise  velocity  gradient,  results  in  a  further  reduction  in  the  skin-friction  drag.  It  is  shown  that 
the  input  power  requirement  is  significantly  less  than  the  power  saved  by  reduced  drag.  Other 
turbulence  characteristics  affected  by  these  controllers  are  also  discussed.  ©  2001  American 
Institute  of  Physics.  [DOI:  10.1063/1.1359420] 


I.  INTRODUCTION 

Much  attention  has  been  paid  to  the  drag  reduction  in 
turbulent  boundary  layers.  Skin  friction  drag  constitutes  ap¬ 
proximately  50%,  90%,  and  100%  of  the  total  drag  on  com¬ 
mercial  aircraft,  underwater  vehicles,  and  pipelines, 
respectively.1  The  decrease  of  skin  friction,  therefore,  entails 
a  substantial  saving  of  operational  cost  for  commercial  air¬ 
craft  and  submarines.  Recent  reviews1'3  summarize  achieve¬ 
ments  and  open  questions  in  boundary  layer  control. 

With  the  notion  that  near-wall  streamwise  vortices  are 
responsible  for  high  skin  friction  in  turbulent  boundary  lay¬ 
ers,  Choi  et  alA  manipulated  the  near-wall  turbulence  by  ap¬ 
plying  various  wall  actuations.  They  achieved  a  20%  skin- 
friction  reduction  in  a  turbulent  channel  flow  by  applying  a 
wall  transpiration  equal  and  opposite  to  the  wall-normal  ve¬ 
locity  component  measured  at  y  +  =  10.  This  control  is  shown 
to  effectively  make  the  streamwise  vortices  weaker.  How¬ 
ever,  it  is  not  easily  implementable  since  it  is  difficult  to 
place  sensors  inside  the  flow  field.  Other  attempts  at  weak¬ 
ening  the  near-wall  streamwise  vortices  have  been  made  by 
imposing  spanwise  oscillation  of  the  wall5  and  using  external 
body  force.6  These  methods,  however,  require  a  large 
amount  of  input  energy.  A  reduction  in  skin  friction  must  be 
accompanied  with  the  required  input  energy  much  less  than 
the  energy  saved  by  the  reduction. 

A  systematic  approach,  not  relying  on  physical  intuition, 
has  been  tried  in  the  past.  A  suboptimal  control,  which  de- 
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termines  the  optimal  control  input  by  minimizing  the  cost 
functional  for  a  short  time  interval,  was  successfully  applied 
to  the  stochastic  Burger’s  equation.7  Bewley  and  Moin8  ex¬ 
tended  the  suboptimal  control  to  a  turbulent  channel  flow. 
This  method,  however,  requires  information  about  the  whole 
flow  field  and  excessive  computation,  so  that  it  is  impossible 
or  at  best  extremely  difficult  to  implement.  It  is  necessary  to 
develop  a  control  scheme  that  utilizes  easily  measurable 
quantities. 

Lee  et  al.9  developed  a  neural  network  control  algorithm 
that  approximates  the  correlation  between  the  wall-shear 
stresses  and  the  wall  actuation  and  then  predicts  the  optimal 
wall  actuation  to  produce  the  minimum  value  of  skin  fric¬ 
tion.  They  also  produced  a  simple  control  scheme  from  this 
neural  network  control,  which  determines  the  actuation  as 
the  sum  of  the  weighted  spanwise  wall-shear  stress, 
dwldy  \w.  Recently,  Koumoutsakos10  reported  a  substantial 
drag  reduction  obtained  by  applying  a  feedback  control 
scheme  based  on  the  measurement  and  manipulation  of  the 
wall  vorticity  flux.  Furthermore,  he  showed  that  the  strength 
of  unsteady  mass  transpiration  actuators  can  be  derived  ex¬ 
plicitly  by  inverting  a  system  of  equations. 

Other  systematic  controls11'17  have  been  developed  by 
exploiting  the  tools  recently  developed  in  the  control 
community.18'21  Joshi  et  a/.11'13  and  Bewley  and  Liu14  de¬ 
veloped  an  integral  feedback  controller,  a  linear  quadratic 
(LQ)  controller,  and  an  H «  controller  (worst-case  controller) 
to  successfully  stabilize  unstable  disturbances  in  transitional 
flow.  In  particular,  Cortelezzi  and  Speyer15  introduced  the 
multi-input-multi-output  (MEMO)  linear-quadratic-Gaussian 
(LQG)/loop-transfer-recovery  (LTR)  synthesis,22  combined 
with  model  reduction  techniques,  for  designing  an  optimal 
linear  feedback  controller.  This  controller  successfully  sup- 
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pressed  near-wall  disturbances,  thus  preventing  a  transition 
in  two-dimensional  laminar  channel  flows.  This  reduced- 
order  controller16  was  applied  to  two-dimensional  nonlinear 
transitional  flows,  illustrating  that  the  controller  designed 
from  the  linear  model  works  remarkably  well  in  nonlinear 
flows. 

Our  purpose  in  the  present  study  is  to  develop  a  realistic 
robust  optimal  controller  that  systematically  determines  the 
wall  actuation,  in  the  form  of  blowing  and  suction  at  the 
wall,  relying  only  on  a  measured  streamwise  velocity  gradi¬ 
ent  to  reduce  skin  friction  in  a  fully  developed  turbulent 
channel  flow.  A  dynamic  representation  of  the  flow  field  is 
required  for  controller  design.  Due  to  the  complexity  and 
nonlinearity  of  the  Navier-Stokes  equations,  it  is  difficult  to 
derive  model-based  controllers.  Therefore,  the  linearized 
Navier-Stokes  equations  for  Poiseuille  flow  are  used  as  an 
approximation  of  the  flow  field  and  form  the  basis  of  system 
modeling.  Several  investigators  (e.g.,  Farrel  and  Ioannou,23 
Kim  and  Lim,24  to  name  a  few)  have  shown  that  linearized 
models  have  a  direct  relevance  to  turbulent  flows.  A  reduced- 
order  controller  has  been  designed  based  on  this  model  and 
applied  to  linear  and  nonlinear  transitional  flows.15-17  En¬ 
couraged  by  these  results,  in  this  paper  we  apply  this  distrib¬ 
uted  two-dimensional  controller  to  a  direct  numerical  simu¬ 
lation  of  turbulent  channel  flow  at  a  low  Reynolds  number. 
We  then  augment  our  two-dimensional  distributed  controller 
by  including  an  ad  hoc  control  scheme  to  attenuate  the  re¬ 
sidual  disturbances  in  the  spanwise  direction. 

In  Sec.  II,  we  derive  the  state-space  equations  from  the 
linearized  two-dimensional  Navier-Stokes  equations.  In  Sec. 
HI,  we  reduce  the  order  of  the  state-space  equations  and 
derive  a  reduced-order  two-dimensional  controller  by  using 
LQG/LTR  synthesis.  In  Sec.  IV,  we  construct  and  apply  the 
distributed  two-dimensional  controller  based  on  the  linear¬ 
ized  Navier-Stokes  equations  to  a  fully  developed  turbulent 
channel  flow  at  Rer-  100,  where  Rer  is  the  Reynolds  num¬ 
ber  based  on  the  wall-shear  velocity,  uT,  and  the  half¬ 
channel  height,  h .  In  Sec.  V,  this  distributed  two-dimensional 
controller  augmented  with  a  simple  ad  hoc  control  scheme  is 
applied  to  the  same  flow.  In  Sec.  VI,  we  present  turbulence 
statistics  associated  with  the  controlled  flows  followed  by 
conclusions  in  Sec.  VII. 

In  this  paper,  we  use  (m,u,w)  to  represent  the  velocity 
components  in  the  streamwise  (jc),  wall-normal  (y),  and 
spanwise  (z)  directions,  respectively. 

il.  THE  STATE-SPACE  EQUATIONS 

One  of  the  goals  in  the  present  study  is  to  reduce  the  size 
of  the  controller.  A  controller  with  a  large  number  of  states  is 
of  no  practical  interest  in  engineering  applications  because  of 
the  amount  of  hardware  and  computer  power  necessary  to 
compute  a  real-time  control  law.  Consequently,  it  is  crucial 
to  reduce  the  order  of  the  controller. 

Figure  1  shows  the  configuration  of  the  turbulent  chan¬ 
nel  flow  equipped  with  the  controller  tested  for  our  study. 
Low-order  controllers  are  usually  preferred  to  high-order  one 
because  of  the  lower  cost  of  hardware  construction  as  well  as 
the  less  computation  time  necessary  to  provide  the  control 


FIG.  1.  Schematic  representation  of  turbulent  channel  flow  equipped  with 
sensors  and  actuators  distributed  in  the  streamwise  direction  in  each  z  plane. 

input.  Hence,  we  slice  the  channel  with  xy  planes  equally 
spaced  in  the  z  direction  in  order  to  reduce  the  order  of  the 
controller.  We  then  construct  the  distributed  two- 
dimensional  controller  by  applying  the  two-dimensional  con¬ 
troller  developed  from  the  linearized  two-dimensional 
Navier-Stokes  equations15  to  each  plane.  It  is  shown16  that 
the  two-dimensional  controller  is  effectively  able  to  reduce 
the  skin-friction  drag  of  the  finite-amplitude  disturbances  in 
a  two-dimensional  channel  flow. 

We  follow  the  same  derivation  of  the  state-space  equa¬ 
tion  as  given  in  Cortelezzi  et  al16  We  give  a  brief  outline 
here  for  completeness;  the  interested  reader  is  referred  to 
Cortelezzi  et  al 16  for  details.  The  wall  transpiration  is  ap¬ 
plied  to  both  top  and  bottom  walls  in  a  fully  developed  tur¬ 
bulent  channel  flow.  For  simplicity,  though,  we  derive  the 
state-space  equations  assuming  that  blowing  and  suction  is 
applied  only  at  the  bottom  wall.  The  application  of  blowing 
and  suction  to  both  walls  is  a  trivial  extension. 

We  consider  two-dimensional  incompressible  Poiseuille 
flow  in  a  periodic  channel  of  streamwise  length,  Lx ,  and 
channel  height,  2  h.  The  undisturbed  velocity  field  has  a 
parabolic  profile  with  centerline  velocity  Uc .  The  linearized 
two-dimensional  Navier-Stokes  equations  can  be  written  in 
terms  of  the  perturbation  streamfunction, 

(d(+tf<?JA(A-t/Vx=Re_1AA^,  (1) 

where  all  variables  are  normalized  with  Uc  and  h  and  Re 
=  Uch!v  is  the  Reynolds  number. 

To  suppress  perturbations  evolving  within  the  bottom 
boundary  layer,  we  apply  blowing  and  suction  at  the  bottom 
wall  (see  Fig.  1).  For  simplicity,  we  assume  that  the  actuators 
are  continuously  distributed.  The  corresponding  boundary 
conditions  are 


Downloaded  14  Sep  2001  to  128.97.82.204.  Redistribution  subject  to  AIP  license  or  copyright,  see  http://ojps.aip.org/phf/phfcr.jsp 


Phys.  Fluids,  Vol.  13,  No.  5,  May  2001 


Application  of  reduced-order  controller  to  turbulent  flows  1323 


0xly=  — 1  ifry  |  y  =  +  j  _y  =  J  (2) 

where  the  control  function  vw  indicates  blowing  and  suction 
at  the  bottom  wall.  We  impose  the  wall  transpiration  of  zero 
net  mass  flux. 

To  detect  the  near-wall  disturbances,  we  measure  the 
gradient  of  the  streamwise  disturbance  velocity  at  given 
point  x=Xj  along  the  bottom  wall  (see  Fig.  1), 

Z(Xi,t)=lpyy  |,—  1.  (3) 

In  other  words,  we  measure  the  first  term  of  the  disturbance 
wall-shear  stress,  Tyx=Rz~\i//yy- ipxx)\y=-X .  The  second 
term  of  the  wall-shear  stress  is  zero  in  the  uncontrolled  case 
and  is  known  in  the  controlled  case. 

We  define  a  performance  index  Jy  or  cost  criterion,  to 
design  a  controller  for  the  LQG  ( H2 )  problem.  Since  we  are 
interested  in  suppressing  the  disturbance  wall-shear  stress, 
ryx ,  we  define 

lim  \S[  (tl'ly+</flx)\y=-\dxdt.  (4) 

/ J - *00  J  t  Jo 

The  integrand  represents  the  cost  of  the  disturbance  wall- 
shear  stress,  ryx ,  being  different  from  zero.  Moreover,  the 
integrand  implicitly  accounts  for  the  cost  of  implementing 
the  control  itself.  There  are  two  reasons  to  minimize  the  cost 
of  the  control.  In  any  engineering  application  the  energy 
available  to  drive  the  controller  is  limited,  and  a  large  control 
action  may  drive  the  system  away  from  the  region  where  the 
linear  model  is  valid. 

By  using  the  same  procedure  described  in  Cortelezzi 
et  ah ,16  Eqs.  (l)-(3)  are  converted  into  the  state-space  equa¬ 
tions: 

dx 

—  =  Ax+Bu,  z=Cx+Du,  (5) 

with  the  initial  condition  x(0)  =  x0,  where  x  is  the  internal 
state  vector,  u  is  the  control  vector,  z  is  the  measurement 
vector.  Matrices  A,B,C  contain  the  dynamics  of  the  two- 
dimensional  plane  Poiseuille  flow,  actuators,  and  sensors,  re¬ 
spectively.  Matrix  D  contains  the  coupling  between  sensors 
and  actuators.  The  cost  criterion,  Eq.  (4),  becomes 

3=  lim  f  (/[  zTz  +  urFrFu]df ,  (6) 

tj—KX>  J  t 

where  the  superscript  T  denotes  a  transposed  quantity.  The 
matrix  F  is  obtained  by  spectrally  decomposing  the  last  term 
in  the  cost  criterion,  Eq.  (4). 

The  advantage  of  the  present  formulation  is  that  the 
whole  problem  decouples  with  respect  to  the  wave  number 
when  Eqs.  (5)  and  (6)  are  transformed  into  Fourier  space  in 
the  streamwise  direction.  All  matrices  in  Eqs.  (5)  and  (6)  are 
block  diagonal,  which  allows  the  above  state-space  system 
into  equivalent  N  state-space  subsystems.25  For  a  given  wave 
number,  a,  the  state-space  subsystem  equations  are 

dxa 

- =  Aaxa+Baua,  za=Caxa+Daiia,  (7) 


with  the  initial  condition  xa(0)  =  xa0.  It  can  be  shown  that 
the  cost  criterion,  Eq.  (6),  also  decouples  with  respect  to  the 
wave  number  (otherwise  the  wave  number  decoupling  is  not 
possible  while  the  system  itself  is  decoupled),  and  we  obtain 
N  performance  indexes.  For  a  given  wave  number,  a,  the 
cost  criterion  is  defined  as 

Ja=  lim  Ja=  lim  [  \zTaza+nTaVTaYa\x a]dt.  (8) 

t  00  tj—>  00  J  t 

Consequently,  the  design  of  a  two-dimensional  controller  for 
the  system,  Eq.  (5),  with  a  specified  cost  criterion,  Eq.  (6), 
has  been  reduced  to  the  independent  design  of  N  single-wave 
number  controllers  for  the  subsystems,  Eq.  (7),  along  with 
Eq.  (8). 


111.  MODEL  REDUCTION  AND  CONTROLLER  DESIGN 

In  this  section  we  derive  a  lower-order  two-dimensional 
controller  in  two  steps.15  First,  we  construct  a  lower-order 
model  of  Eq.  (7),  and  subsequently,  design  a  single-wave 
number  controller  for  the  reduced-order  model.  To  obtain  a 
lower-order  model,  we  transform  Eq.  (7)  into  a  Jordan  ca¬ 
nonical  form.  The  matrices  Aa,Ba,C„,Da  that  describe  the 
dynamics  of  the  reduced-order  model  are  obtained  from  the 
matrices,  Aa,Ba,Ca,Da  in  the  Jordan  canonical  form  by 
retaining  rows  and  columns  corresponding  to  equally  well 
controllable  and  observable  states.  The  overcaret  denotes  the 
quantities  associated  with  the  reduced-order  model. 

Although  a  rigorous  mathematical  framework  for  the  de¬ 
sign  of  disturbance  attenuation  (Tf*)  linear  controllers  is 
provided  by  the  control  synthesis  theory,18,19  for  this  study 
LQG(7f2)  synthesis  is  preferred.  A  brief  review  will  be 
given  in  a  self-contained  manner  to  provide  the  necessary 
governing  equations  for  closed-loop  stability  analysis.20 

The  LQG  problem  for  each  wave  number  a  is  formu¬ 
lated  as  a  stochastic  optimal  control  problem  described  by 
equations 

x*= Aaxa+B„ua+  r„wa ,  (9) 

z«=C*xa+Daua+va,  (10) 

where  Fa  is  an  input  matrix,  wa  and  \a  are  both  white  noise 
processes  with  zero  means  and  autocorrelation  functions, 

£[wa(f)w£(r)]  =  W„<S(f-T), 

X  (11) 

£[va(f)v£(T)]  =  V„<5(f-T), 

where  £[  •  ]  is  the  expectation  operator  averaging  over  all 
underlying  random  variables  and  S(t-r)  is  the  delta  func¬ 
tion.  Note  that  Wa  and  Va ,  the  power  spectral  densities,  will 
be  chosen  here  as  design  parameters  to  enhance  system  per¬ 
formance.  An  additional  comment  on  the  controller  design 
process  will  be  given  at  the  end  of  this  section. 

The  LQG  controller  is  determined  by  finding  the  control 
action  u a(Z,),  where  Zf  =  {z(r);0^  r^f}  is  the  measure¬ 
ment  history,  which  minimizes  the  cost  criterion, 
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1 

Ja—  lim  - 

tj- +00 

t  /(xlQaxa+  2x^Naua+  uTaRaua)d^j ,  (12) 

subject  to  the  stochastic  dynamics  system  model  equations 
Eqs.  (10)— (1 1).  Note  that,  from  Eqs.  (7)-(8),  Q*=C£Ca, 
Na=C^Da,  and  Ra=D^Da+F^Fa .  The  division  by  ( tf 
—  t)  ensures  that  the  cost  criterion  remains  finite  in  the  pres¬ 
ence  of  uncertainties  in  the  infinite-time  problem  {tf- *o°). 
Note  that  Eq.  (12)  can  include  Eq.  (8),  where 

1 

Ja=  lim  - E[Ja ],  (13) 

tj—> »  tf  t 

and  the  limit  in  Eq.  (8)  is  explicitly  denoted  in  Eq.  (13).  Note 
that  even  though  the  time  interval  is  infinite,  the  time  re¬ 
sponse  is  still  measured  by  the  eigenvalues  of  the  closed- 
loop  system.  We  consider  the  infinite-time  problem  with 
time-invariant  dynamics  because  the  controller  gains  become 
constants. 

By  nesting  the  conditional  expectation  with  respect  to  Zt 
within  the  unconditional  expectation  of  Eq.  (13),  i.e., 
E[Jai  —  E[E[Ja/Zt]]9  where  E[-IZt]  denotes  the  expecta¬ 
tion  ( • )  conditioned  on  Zt ,  the  cost  criterion  can  be  written 
as 


/  =  lim 


tj— if 


tf-t 


x£(  //^«Q“^+2^N“u“+u«R“u“+tr(p“^rfT)  > 


(14) 

where  xa—E[xa/Zt]  is  the  conditional  mean  estimate  of  the 
state  x  and  Pa  is  the  conditional  error  variance.  This  cost 
criterion  is  now  minimized  subject  to  the  estimation  equa¬ 
tions  discussed  below.  Note  that  Pa  does  not  depend  on  the 
control  [see  Eq.  (18)  below]  and,  therefore,  does  not  enter 
into  the  optimization  process. 

The  solution  to  the  regulator  problem20  is  a  compensator 
composed  of  a  state  reconstruction  process,  known  here  as  a 
filter  (in  the  no-noise  case  it  is  known  as  an  observer)  in 
cascade  with  a  controller  (see  Fig.  1,  where  £,  is  the  estima¬ 
tor  and  Cj  is  the  controller).  The  state  estimate  (conditional 
mean)  xa  is  governed  by  the  so-called  Kalman  filter  as 


xa  =  A  axa + B  aua  +LavaJ 
va  =  za  - za  =  Ca(xa-xa)  +  \a . 


(15) 


If  the  reduced-order  system  were  the  actual  system,  then  va 
in  Eq.  (15)  is  correct.  When  the  actual  system  is  considered 
and  the  filter  is  implemented  based  on  the  reduced-state 
space,  z  rather  than  z  is  the  measurement  and  the  filter  re¬ 
sidual  becomes 


,=  za-Caxtt-D„ua 


(16) 


The  Kalman  gain  matrix  La ,  constructed  to  trade  the  accu¬ 
racy  of  the  new  measurements  against  the  accuracy  of  the 
state  propagated  from  the  system  dynamics,  is  given  by 

L«=P^Va"-  (17) 

where  Pa  is  the  error  variance  in  the  statistical  problem. 

In  the  infinite-time  stationary  formulation,  the  error  Fa  is 
the  solution  to  the  algebraic  Riccati  equation  (ARE), 

Aapa+ paA^+ rffwar^-  p„cX 1  capa= o.  (is) 

If  the  system  is  (Aa,Ca)  observable  and  (Aa,Ba)  control¬ 
lable,  then  Pa  is  positive  definite.  Under  these  assumptions, 
it  can  be  shown  that  the  difference  between  the  internal  state 
xa  and  the  estimate  state  xa ,  i.e.,  the  error, 

ea=xa-xa ,  (19) 

goes  to  zero  as  time  goes  to  infinity.  In  other  words,  the 

evolution  equation, 

ea= A^+  \jav+  Taw,  (20) 

is  stable,  i.e.,  all  the  eigenvalues  of  the  matrix, 

A/=Aa— LaCa ,  (21) 

have  a  negative  real  part. 

Minimizing  the  infinite-time  cost  function  /,  Eq.  (14) 
subject  to  Eq.  (15)  yields  the  following  control  law: 

u*=-Kaxa,  (22) 

where 


K^R^BX+NJ,  (23) 

and  Sa  is  the  solution  of  the  algebraic  Riccati  equation 
(ARE), 

A„S„+  S„ &Ta+  Qa~(SaBa+ N„)R“ 1  (BX+ N«)  =  0.  (24) 

It  should  be  remarked  that  the  control  gain  matrix  Ka  is 
determined  from  functions  only  of  the  known  dynamics  co¬ 
efficients  (Aa,Ba)  and  the  weighting  in  the  cost  criterion 
(Qa,RJ,  and  not  the  statistic  of  the  input  (Va,  Wa).  Con¬ 
sequently,  Ka  is  determined  from  a  performance  index  as 
Eq.  (12),  independent  of  the  stochastic  inputs.  If  (Aa,Ba)  is 
controllable  and  (A^jQ^2)  observable,  then  the  loop  coeffi¬ 
cient  matrix, 

Ac=As-KaBa,  (25) 

is  stable  and  Sa  is  positive  definite.  The  controllable  and 
observable  conditions  can  be  weakened  to  stabilizable  and 
detectable.21 

When  we  combine  the  estimator  and  the  regulator  to¬ 
gether,  the  dynamic  system  composed  of  the  controlled  pro¬ 
cess  and  filter  becomes 


0 

Ar 


M+(L<*v“+r“w“ 

xj  \  taya 


(26) 


Note  that  any  choice  of  two  between  e,  xa ,  and  xa  produce 
the  same  dynamics  because  they  are  algebraically  related  by 
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Eq.  (19).  Under  the  above  controllability  and  observability 
assumptions,  Ay  and  Ac  have  only  stable  eigenvalues  if  op¬ 
timal  gains  La  and  Ka  of  Eqs.  (17)  and  (23)  are  used.  If  the 
actual  linear  system  is  used,  then  xa  and  the  reduced-order 
state  estimate  xa  are  used  to  form  the  closed-loop  dynamic 
system  rather  than  that  given  in  Eq.  (26).  The  eigenvalues  of 
the  dynamical  matrix  now  dictate  the  system  stability  and 
will  differ  from  the  ideal  case  of  Eq.  (26). 

The  parameters  used  in  our  LQG  design  are  now  ad¬ 
dressed.  Since  the  power  spectral  density  is  not  known,  for 
simplicity  of  the  design  we  consider  Va  and  Wa  to  be  of  the 
form  \a= 01  and  Wa= pi,  where  fi  and  p  are  scalar  and  I  is 
an  identity  matrix.  Only  the  ratio  of  /3  and  p  is  important. 
Furthermore,  by  choosing  ra=Ba,  loop- transfer  recovery 
(LTR)  of  the  LQG  controller  to  full-state  feedback18  guaran¬ 
tees  that  robust  performance  occurs  when  the  process  noise 
power  spectral  density  goes  to  infinity,  i.e.,  p— >oo,  provided 
there  exists  no  nonminimal-phase  zero  in  the  plant.  In  our 
case,  there  are  no  nonminimal-phase  zeros  and  robust  perfor¬ 
mance  means  approximately  obtaining  60°  of  phase  margin 
and  6  db  of  the  gain  margin.  Note  that  the  choice  of  Ta 
=  Ba  implies  that  the  noise  is  generated  along  the  wall  as  is 
the  control  and  could  be  interpreted  as  due  to  wall  roughness. 
Furthermore,  the  values  of  p  and  p  were  determined  by  tun¬ 
ing  the  controller  in  the  presence  of  turbulent  flow.  The  de¬ 
gree  of  loop  transfer  recovery  varied  from  controller  to  con¬ 
troller. 

As  described  above  using  LQG/LTR  assumes  that  the 
uncertainty  is  at  the  wall  and  effects  the  dynamics  in  the 
same  way  as  the  control.  Since  the  system  has  the  same 
controllability  with  respect  to  both  the  control  and  distur¬ 
bances,  state-space  reduction  for  controller  design  was 
straightforward.  This  is  in  contrast  to  Hoo  control  used  by 
Bewley  and  Liu,14  where  uncertainty  is  assumed  uniformly 
throughout  the  channel.  Since  controllability  of  the  distur¬ 
bances  is  different  from  that  of  the  control,  model  reduction 
may  not  be  straightforward.  Furthermore,  robustness  in  terms 
of  traditional  measures  of  the  gain  and  phase  margin  in  con¬ 
trol  engineering  are  also  obtained  by  using  LQG/LTR.  For 
these  reasons  LQG/LTR  is  used  for  the  present  study  instead 
of  the  unstructured  uncertainty  controllers. 

Figure  1  links  the  mathematical  formulation  to  its  com¬ 
putational  implementation  by  summarizing  in  a  block  dia¬ 
gram  the  control  strategy  described  above.  The  two- 
dimensional  distributed  controller  can  be  programmed  in  a 
computer  routine  whose  input  is  a  matrix  containing  the  gra¬ 
dients  of  the  streamwise  velocity  component  and  whose  out¬ 
put  is  a  matrix  containing  the  blowing  and  suction  at  the 
wall.  Each  column  of  the  measurement  matrices  contains  the 
gradients  of  the  streamwise  velocity  component  along  the 
wall  at  a  given  spanwise  location.  Each  column  is  processed 
in  parallel  by  a  fast  Fourier  transform  (FFT)  and  converted 
into  za’s.  Each  single-wave  number  controller,  Eqs.  (9)- 
(10),  is  integrated  in  time  by,  for  example,  a  third-order  low- 
storage  Runge-Kutta  scheme.  The  ua’s  are  computed  in  par¬ 
allel.  An  inverse  FFT  converts  ua’s  into  the  columns  of  the 
matrix  containing  the  blowing  and  suction  at  the  wall  along 
the  streamwise  direction.  This  routine  can  be  embedded  in 
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any  Navier-Stokes  solver  able  to  handle  time-dependent 
boundary  conditions  for  the  control  of  three-dimensional 
channel  flows. 

Figure  1  also  provides  the  basic  architecture  for  the  po¬ 
tential  implementation  of  the  present  distributed  two- 
dimensional  controller  in  practical  engineering  applications. 
For  instance,  the  gradients  of  the  streamwise  velocity  com¬ 
ponent  can  be  measured  by  microelectromechanical-systems 
(MEMS)  hot-film  sensors.26  For  each  xy  plane,  analog-to- 
digital  converters  (A/D)  and  digital  signal  processors  (DSP) 
convert  the  measured  gradients  into  ztt’s.  Each  single- wave 
number  controller,  Eqs.  (9)-(10),  is  replaced  by  a  micropro¬ 
cessor,  and  parallel  computation  produces  ua’s.  A  DSP  and  a 
digital-to-analog  converter  (D/A)  produce  the  actuating  sig¬ 
nal  in  each  xy  plane.  A  variety  of  actuators,  such  as  synthetic 
jets,  microbubble  actuators,  and  thermal  actuators,  can 
mimic  small-amplitude  blowing  and  suction  at  the  wall.26 

Although  the  structure  of  this  compensator  is  simplified 
by  the  parallel  computation  (for  all  spanwise  directions),  it 
does  require  processing  of  all  the  sensor  measurements  (for 
all  streamwise  directions).  The  controller  is  essentially  cen¬ 
tralized  because  all  information  is  used  and  the  actuators  are 
activated  spatially  over  the  assumed  channel.  Controllers 
based  explicitly  on  the  spatial  distribution  of  the  control, 
suggested  by  Bamieh  et  al}1  show  that  there  is  a  spatial 
decay  rate.  Our  controller  can  be  constructed  to  represent  a 
discrete  form  of  their  controller  and  given  the  spatial  decay 
rate  for  our  configuration,  i.e.,  the  size  of  the  channel  could 
be  chosen  consistent  with  that  decay  rate.  Nevertheless,  our 
representation  allows  a  significant  decrease  in  on-line  com¬ 
putation  by  identifying  the  Fourier  modes  and  the  number  of 
states  associated  with  those  modes  that  best  reduce  turbu¬ 
lence  as  discussed  in  the  next  section. 

IV.  PERFORMANCE  OF  A  TWO-DIMENSIONAL 
CONTROLLER 

For  the  purpose  of  testing  the  performance  of  a  control¬ 
ler,  we  performed  direct  numerical  simulations  of  a  turbulent 
channel  flow  at  Rer=  100.  A  spectral  code  was  used  with  a 
computational  domain  of  (47r,2,47r/3)  and  a  grid  resolution 
of  (32,65,32)  in  the  (x,y,z)  directions.  The  numerical  tech¬ 
nique  used  in  this  study  is  essentially  the  same  as  that  of  Kim 
et  al2%  except  that  the  time  advancement  for  the  nonlinear 
terms  is  a  third-order  Runge-Kutta  (RK3)  method.  The 
second-order  accurate  Crank-Nicolson  (CN)  method  is  used 
for  the  linear  terms. 

We  designed  a  distributed  two-dimensional  controller  in 
two  steps.  First,  we  designed  reduced-order  controllers  for 
two-dimensional  Poiseuille  flow  in  a  periodic  channel  of 
streamwise  length  Lx-Att  dX  Re=5000,  which  has  the  same 
mean  wall-shear  stress  as  turbulent  channel  flow  at  Rer 
=  100.  Subsequently,  we  fine-tuned  the  single-wave  number 
reduced-order  controllers  in  order  to  minimize  the  magnitude 
of  the  Fourier  coefficients  of  the  wall-shear  stresses  in  tur¬ 
bulent  channel  flow  at  ReT=  100.  We  used  N=  32  and  M 
=  60  in  this  linear  model  flow.  Controllers  operate  at  both 
top  and  bottom  walls  in  parallel.  If  the  two-dimensional  con¬ 
trollers  without  model  reduction  were  applied  at  each  z 
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FIG.  2.  Time  history  of  the  drag  for  the  controlled  and  uncontrolled  flows: 
- ,  controlled  flow; - ,  uncontrolled  flow. 

plane,  then  the  order  of  the  ensemble  of  controllers  would  be 
64X3904=249  856.  Using  the  model  reduction  technique 
previously  described,  we  designed  eight  single-wave  number 
controllers  of  order  12,  corresponding  to  the  eight  lowest 
wave  numbers.  Since  we  use  the  eight  lowest  single- wave 
number  controllers  in  our  simulation,  the  combined  order  of 
the  controllers  is  64X96=6144.  It  represents  a  state-space 
reduction  of  about  97.5%,  with  respect  to  the  full-order  sys¬ 
tem. 

Figure  2  shows  the  time  history  of  the  drag  in  the  un¬ 
controlled  and  controlled  flows.  Drag  is  measured  by  the 
mean  value  of  the  wall-shear  stresses  averaged  over  each  top 
and  bottom  wall.  This  two-dimensional  control  yields  about 
a  10%  drag  reduction.  Choi  et  alA  reported  that  the  in-phase 
u  control  measured  at  y  +  =  10  also  gives  a  10%  drag  reduc¬ 
tion.  This  in-phase  streamwise  velocity  at  the  wall  causes  a 
similar  effect,  du'ldy\w— 0,  which  is  the  to-be-minimized 
target  of  our  cost  criterion  in  our  two-dimensional  controller. 
Note  that  this  observed  drag  reduction  is  a  byproduct  since 
our  controller  is  designed  to  suppress  the  fluctuations  of  the 
streamwise  wall-shear  stress,  not  the  mean  wall-shear  stress. 
Note  also  the  sudden  drop  in  the  drag  as  soon  as  the  control¬ 
ler  is  switched  on  at  t= 25.  This  transient  phenomena  is  also 
observed  in  other  studies.8,9 

Figure  3  compares  the  magnitude  of  Fourier  coefficients 
of  the  wall-shear  stresses  in  the  controlled  and  uncontrolled 
flows.  The  wall-shear  stresses  are  measured  at  the  bottom 
wall  at  a  given  spanwise  location.  Figures  3(a)  and  3(b)  show 
the  comparisons  corresponding  to  wave  numbers  kx=  0.5  and 
kx— 1.0,  respectively.  Both  figures  show  an  order-of- 
magnitude  reduction  between  the  controlled  and  uncon¬ 
trolled  cases.  The  magnitude  of  the  Fourier  coefficients  of 
wall-shear  stress  decreases  very  quickly  as  soon  as  the  con¬ 
troller  is  activated  at  r=25.  These  results  indicate  that  our 
distributed  two-dimensional  linear  reduced-order  controller 
suppresses  disturbance  wall-shear  stress  remarkably  well, 
even  in  a  fully  developed  turbulent  flow.  The  high  wave 
number  components  of  the  wall-shear  stress  in  Fig.  3(c)  do 
not  show  any  reduction  since  only  the  lowest  eight  single¬ 
wave  number  controllers  (up  to  kx— 4.0)  are  used  in  the  con¬ 
trol  of  flow.  Examinations  of  other  spanwise  locations  show 
similar  results. 

Contours  of  the  disturbance  wall-shear  stresses  at  the 


FIG.  3.  Time  history  of  the  magnitude  of  the  Fourier  coefficients  of  the 
wall-shear  stresses  measured  at  the  bottom  wall  at  a  given  spanwise  location 

for  the  controlled  and  uncontrolled  flows: - ,  uncontrolled  flow; - , 

controlled  flow,  (a)  £*=0.5,  (b)  kx~  1.0,  and  (c)  £*=6.0. 


bottom  wall  in  the  controlled  and  uncontrolled  flows  at  t 
=  30  are  shown  in  Fig.  4.  Contours  for  the  uncontrolled  flow 
show  the  usual  elongated  regions  of  low-  and  high-shear 
stress.  Note  that  contours  for  the  controlled  flow  show  the 
dramatic  effect  of  the  distributed  two-dimensional  controller. 
The  long  streaky  wall-shear  stress  region  spans  almost  the 
entire  streamwise  direction,  indicating  that  the  low  wave 
number  components  (except  the  zero  wave  number  that  we 
do  not  control)  are  completely  suppressed,  which  is  consis¬ 
tent  with  Fig.  3.  The  remaining  spanwise  variations,  i.e.,  the 
alternating  regions  of  high-  and  low-shear  stress,  are  due  to 


X 


FIG.  4.  Contours  of  disturbance  wall-shear  stresses  at  the  bottom  wall  at  t 
=  30:  (a)  uncontrolled  flow;  (b)  2-D-controlled  flow.  Negative  contours  are 
dashed. 
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the  fact  that  the  two-dimensional  controllers  distributed 
along  the  streamwise  direction  are  operated  independently 
from  one  z  plane  to  another. 

The  above  results  demonstrate  that  our  distributed  two- 
dimensional  controller  designed  from  the  linear  model  works 
remarkably  well  in  suppressing  near-wall  disturbances  in  the 
fully  developed  turbulent  flow.  The  reduction  of  fluctuating 
wall-shear  stress  led  to  drag  reduction.  However,  this  distrib¬ 
uted  two-dimensional  controller  has  a  limited  impact  on  the 
total  drag  since  it  cannot  control  the  spanwise  variation  of 
the  wall-shear  stress.  In  the  next  section  an  augmentation  to 
the  distributed  two-dimensional  controller  is  presented  and 
implemented. 

V.  AUGMENTED  THREE-DIMENSIONAL  CONTROLLER 

In  the  previous  section,  successful  control  of  fully  devel¬ 
oped  turbulent  channel  flow  has  been  obtained  by  applying  a 
distributed  two-dimensional  controller.  However,  it  has  been 
observed  that  this  controller  does  not  account  for  the  span- 
wise  variations  of  fluid  motion.  An  augmentation  to  the  dis¬ 
tributed  two-dimensional  controller  that  accommodates  the 
three-dimensional  characteristics  of  a  fully  developed  turbu¬ 
lent  flow  is  developed  in  this  section. 

A  simple  ad  hoc  control  augmentation  scheme  is  intro¬ 
duced  in  an  attempt  to  capture  the  remaining  spanwise  varia¬ 
tions  of  the  controlled  flow.  This  additional  control,  which 
generates  blowing/suction  to  attenuate  the  spanwise  variation 
of  the  wall-shear  stress,  is  given  as  follows: 

Idu  {x'z)  du  x\ 

e7> 

where  duldy\^x'z)  and  dufdy\xw  are  the  streamwise  velocity 
gradients  averaged  over  the  x z  plane  and  the  x  direction, 
respectively,  and  C  is  a  constant  to  be  adjusted  for  the  best 
performance.  The  subscript  ad  indicates  the  ad  hoc  control, 
and  uad  ts  a  function  of  only  z.  Therefore,  the  new  control 
input  is  defined  by 

vw(xyz)  =  vaLd+  u2-d,  (28) 

where  u2-d  is  the  actuation  velocity  generated  by  the  distrib¬ 
uted  two-dimensional  controller  used  in  the  previous  section. 

Using  the  distributed  two-dimensional  controller  aug¬ 
mented  with  this  ad  hoc  control  scheme,  the  control  of  the 
fully  developed  turbulent  flow  with  Rer=  100  increased  drag 
reduction  to  about  17%,  as  shown  in  Fig.  5.  As  before,  the 
turbulent  flow  is  left  free  to  evolve  without  any  wall  actua¬ 
tion  until  /  =  25.  As  soon  as  the  controller  is  activated  at  t 
=  25,  the  drag  drops  sharply  within  a  very  small  time  period. 
The  constant,  C,  in  Eq.  (27)  is  adjusted  such  that  the  root- 
mean-square  (rms)  value  of  the  actuation  is  maintained  at 
0.1  wr,  where  uT  is  the  wall-shear  velocity  for  the  uncon¬ 
trolled  flow.  We  have  found  empirically  that  C  between 
0.05wr  and  0.2wr  gives  a  similar  performance.  An  introduc¬ 
tion  of  this  simple  control  augmentation  enhances  the  drag 
reduction,  indicating  that  more  sophisticated  controllers  that 
best  take  into  account  the  three-dimensionality  of  turbulent 
flow  may  produce  even  more  efficient  suppression  of  skin- 
friction  drag. 
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FIG.  5.  Time  history  of  the  drag  for  the  controlled  and  uncontrolled  flows: 

- ,  uncontrolled  flow; - ,  2-D-controlIed  flow;  •  ■  *,  ad  /ioc-controlled 

flow. 


Figure  6  presents  the  comparison  of  contours  of  the  dis¬ 
turbance  wall-shear  stresses  at  the  bottom  wall  between  the 
ad  hoc  controlled  flow  and  the  uncontrolled  flow  at  /  =  30. 
Compared  to  Fig.  4,  additional  effort  in  the  spanwise  direc¬ 
tion,  v  acj ,  removes  the  pronounced  peak- valley  variation  of 
the  wall-shear  stress  that  is  observed  in  the  controlled  flow 
with  the  distributed  two-dimensional  controllers  [see  Fig. 
4(b)].  Note  that  the  high  wave  number  components  of  the 
wall-shear  stress  are  persistently  sustained  because  of  the 
lowest  eight  single-wave  number  controllers  adopted  in  the 
control  of  flow. 


VI.  TURBULENCE  STATISTICS 

Some  turbulence  statistics  of  the  flow  field  associated 
with  the  two  controllers  applied  in  this  paper  were  examined 
to  investigate  the  effect  of  the  controllers  on  turbulence.  All 
statistical  quantities  were  averaged  over  a  sufficiently  long 
interval  of  time  as  well  as  over  the  planes  parallel  to  the  wall. 
For  simplicity,  the  flows  controlled  by  the  distributed  two- 
dimensional  controller  only  and  the  distributed  two- 
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FIG.  6.  Contours  of  disturbance  wall-shear  stresses  at  the  bottom  wall  at  t 
=  30:  (a)  uncontrolled  flow;  (b)  ad  ftoc*controlled  flow.  Negative  contours 
are  dashed. 
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FIG.  7.  Mean- velocity  profiles:  ad  /wc-controlled  flow; - ,  2-D- 

controlled  flow; - ,  uncontrolled  flow. 

dimensional  controller  augmented  with  the  ad  hoc  control 
scheme  are  called  ‘^-D-controlled*  ’  and  “ ad  hoc- 
controlled”  flows,  respectively. 

The  mean  velocity  profiles  normalized  by  the  actual 
wall-shear  velocities  are  shown  in  Fig.  7  for  three  different 
channel  flows.  These  profiles  show  the  same  trend  shown  in 
the  Choi  et  alA  drag-reduced  flow:  the  slope  of  the  log  law 
for  controlled  flows  remains  the  same  while  the  mean  veloc¬ 
ity  itself  is  shifted  upward  in  the  log-law  region. 

The  root-mean-square  (rms)  values  of  turbulent  velocity 
fluctuations  are  shown  in  Fig.  8  and  compared  to  those  of  the 


3 


FIG.  8.  Root-mean-square  values  of  turbulent  velocity  fluctuations  normal¬ 
ized  by  the  wall-shear  velocity,  uT  for  the  uncontrolled  flow: - ,  un¬ 
controlled  flow; - ,  2-D-controlled  flow;  •  •  *,  ad  hoc -controlled  flow. 


FIG.  9.  Root-mean-square  values  of  vorticity  fluctuations  normalized  by  the 

wall-shear  velocity  in  wall  coordinates: - ,  uncontrolled  flow; - , 

2-D-controlled  flow;  •  *  \  ad  /zoc-controlled  flow. 


uncontrolled  flow.  Note  that  all  quantities  in  this  figure  are 
normalized  by  the  wall-shear  velocity  of  the  uncontrolled 
flow.  The  controllers  reduce  the  value  of  turbulent  intensity 
significantly  throughout  the  channel,  especially  for  the  wall- 
normal  and  spanwise  components.  The  reduction  of  these 
quantities  in  the  ad  /zoc-controlled  flow  is  greater  than  that  in 
the  2-D-controlled  flow.  The  increase  in  very  near  the 
wail  is  due  to  the  control  input.  A  similar  feature  is  also 
observed  by  Choi  et  al.4  and  Lee  et  al.9  Both  controllers 
mitigate  the  rms  of  spanwise  velocity  fluctuation  throughout 
the  channel  compared  to  that  in  uncontrolled  flow.  However, 
the  introduction  of  uad  in  Eq.  (27)  causes  this  value  to  in¬ 
crease  very  close  to  the  wall,  which  also  leads  to  an  increase 
in  the  streamwise  vorticity  at  the  wall. 

Root-mean-square  values  of  vorticity  fluctuations  for  the 
controlled  flows  are  compared  with  those  for  the  uncon¬ 
trolled  flow  in  Fig.  9.  All  components  of  vorticity  fluctua¬ 
tions  are  significantly  reduced  throughout  the  channel.  Very 
close  to  the  wall,  however,  the  increase  of  streamwise  vor¬ 
ticity  in  the  ad  hoc- controlled  flow  is  due  to  the  streamwise 
vorticity  built  at  the  wall  by  the  ad  hoc  controller.  The  high 
streamwise  vorticity  at  the  wall  slows  the  sweeping  motion 
of  high-momentum  fluid  induced  by  the  streamwise  vorticity 
away  from  the  wall,  thus  resulting  in  a  significant  reduction 
in  skin  friction.  A  similar  feature  is  also  observed  in  Lee 
et  al9  Note  that  the  streamwise  vorticity  at  the  wall  for  the 
2-D-controlled  flow,  however,  is  less  than  that  for  the  uncon- 
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FIG.  10.  A  comparison  of  streamwise  vorticity  contours  in  a  yz  plane  be¬ 
tween  controlled  and  uncontrolled  flows:  (a)  uncontrolled  flow;  (b)  2-D- 
controlled  flow;  (c)  ad  /ioc-controlled  flow.  Negative  contours  are  dashed. 

trolled  flow.  The  reduction  of  <oz  is  a  direct  consequence  of 
the  controller,  which  was  designed  to  reduce  du'ldy\w .  The 
reduction  of  (oy  also  indicates  that  our  controllers  weaken  the 
strength  of  near-wall  streaks.  This  also  decreases  the  streak 
instability,  which  is  shown  to  be  responsible  for  regenerating 
the  near-wall  streamwise  vortices.29,30 

Figure  10  compares  the  streamwise  vorticity  fields  in  the 
uncontrolled  and  controlled  flows.  The  strength  of  the  near¬ 
wall  streamwise  vorticity  for  the  controlled  flows  are  greatly 
attenuated  due  to  the  wall  transpiration  produced  by  the  con¬ 
trollers.  It  is  discernible  that  the  ad  hoc  controller  diminishes 
the  streamwise  vorticity  substantially  more.  The  reduction  of 
the  strength  of  the  streamwise  vorticity  has  also  been  ob¬ 
served  by  Lee  et  ai9  While  Lee  et  al9  suppressed  the 
streamwise  vorticity  field  with  the  physical  understanding 
that  the  control  based  on  the  weighted  sum  of  dw!dy\w  can 
prevent  the  physical  eruption  at  the  wall,  the  present  control¬ 
lers  attenuate  the  streamwise  vorticity  strength  by  minimiz¬ 
ing  the  streamwise  disturbance  wall-shear  stress  systemati¬ 
cally.  The  present  results  further  support  the  notion  that  a 
successful  attenuation  of  the  near- wall  streamwise  vortices 
results  in  a  significant  reduction  in  skin-friction  drag.4 

VII.  CONCLUSIONS 

A  reduced-order  linear  feedback  control  based  on  a  dis¬ 
tributed  two-dimensional  controller  design  is  applied  to  a 
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turbulent  channel  flow.  A  controller  based  on  a  reduced 
model  of  the  linearized  Navier- Stokes  equations  for  a  lami¬ 
nar  Poiseuille  flow  was  designed  by  using  LQG  (H2)/LTR 
synthesis.  This  controller  was  implemented  using  input  mea¬ 
surements  that  are  the  gradients  of  the  streamwise  distur¬ 
bance  velocity  and  output  controls  that  are  the  blowing  and 
suction  at  the  wall. 

First,  we  applied  the  distributed  two-dimensional  con¬ 
troller  to  both  walls  of  a  turbulent  channel  flow  at  ReT 
=  100.  Eight  single-wave  number  controllers  corresponding 
to  eight  lowest  wave  numbers,  reducing  the  order  of  the  con¬ 
troller  about  2.5%  of  the  order  of  the  full  size  system,  are 
applied  to  attain  a  skin-friction  reduction  of  10%  with  re¬ 
spect  to  the  uncontrolled  turbulent  flow.  Next,  a  simple  ad 
hoc  augmented  control  scheme  of  the  distributed  two- 
dimensional  controller  is  introduced  to  capture  the  three- 
dimensionality  of  turbulent  flow.  The  control  of  fully  devel¬ 
oped  turbulent  flow  by  the  distributed  two-dimensional 
controller  augmented  by  the  ad  hoc  control  scheme  produces 
a  17%  reduction  in  skin-friction  drag.  Motivated  by  this  re¬ 
sult,  we  are  currently  developing  controllers  to  more  effi¬ 
ciently  account  for  the  three-dimensionality  of  turbulent 
flow. 

It  should  be  noted  that  the  present  controller,  which  is 
based  on  a  reduced-order  linear  system,  has  achieved  its  de¬ 
sign  objective,  i.e.,  minimization  of  the  wall-shear  stress  dis¬ 
turbances,  quite  remarkably  when  applied  to  the  nonlinear 
flow.  It  was  anticipated  that  the  reduction  of  disturbances 
would  also  lead  to  a  substantial  reduction  of  the  mean  wall- 
shear  stress.  Unfortunately,  this  turned  out  not  to  be  the  case, 
suggesting  that  some  other  cost  functions  should  be  ex¬ 
plored.  By  comparing  with  our  previous  results,9,31  it  was 
found  that  the  present  controller  is  not  as  effective  in  dimin¬ 
ishing  the  strength  of  the  streamwise  vortices  in  the  buffer 
layer,  which  was  the  primary  target  for  other  controllers,  but 
achieved  its  design  goal  by  mainly  affecting  the  region  very 
close  to  the  wall.  In  this  regard,  minimization  of  the  total 
disturbance  energy  in  the  flow  field32  or  minimization  of  the 
linear  coupling  term24  appears  to  be  a  good  candidate  to  be 
explored.  Whether  either  of  these  cost  criterion  is  indeed 
controllable  in  nonlinear  flows,  however,  remains  to  be  in¬ 
vestigated. 

This  study  is  carried  out  at  low  Reynolds  number. 
Whether  our  controller,  based  on  the  reduced-order  linear 
model,  would  work  in  other  turbulent  flows,  should  be  drawn 
from  real  experiments  or  simulations  at  high  Reynolds  num¬ 
ber.  However,  we  expect  that  it  should  work  equally  well  for 
high  Reynolds  number  flow  since  our  controller,  derived 
from  LQG/LTR  synthesis,  recovers  the  robustness  of  LQR, 
whose  characteristics  have  been  partially  tested  over  the  dif¬ 
ferent  Reynolds  number  flows.33 

The  statistics  of  controlled  and  uncontrolled  flows  are 
compared.  The  mean  velocity  profile  is  shifted  upward  in  the 
log  region,  a  typical  characteristic  of  drag-reduced  flow.  Ve¬ 
locity  and  vorticity  fluctuations  as  well  as  Reynolds  shear 
stress  (not  shown)  are  significantly  reduced  due  to  the 
blowing/suction  generated  by  the  controller.  However,  a  ma¬ 
jor  change  is  confined  to  the  wall  region.  Instantaneous  flow 
fields  show  that  the  distributed  two-dimensional  controller 
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attenuates  and  modifies  the  streaky  structure  of  the  boundary 
layer.  Streaks  are  observed  to  span  the  entire  streamwise 
direction  with  velocity  variations  in  the  spanwise  direction. 
These  variations  are  substantially  reduced  by  the  augmented 
controller. 

The  three-dimensional  aspect  of  the  distributed  two- 
dimensional  controller  by  the  augmentation  of  the  ad  hoc 
control  further  reduced  the  skin-friction  drag.  This  three- 
dimensional  controller  produces  secondary  streamwise  vor- 
ticity  at  the  wall,  which  slows  the  sweeping  motions  of  high- 
momentum  fluid  induced  by  the  streamwise  vorticity  away 
from  the  wall.  This  induced  retarding  of  the  primary  stream- 
wise  vorticity  leads  to  additional  drag  reduction,  which  was 
also  observed  in  Choi  et  alA 

Regarding  the  scaling  factor  C  in  Eq.  (27),  we  found  an 
optimal  value  of  C  that  yields  the  blowing/suction  of  0.1m  T. 
With  this  optimal  C,  the  augmented  controller  generates  wall 
transpiration  with  a  rms  value  of  about  0.12mt.  The  required 
power  input  per  unit  area  to  the  system,  pwvw+0.5pvl 
**0.1  pu\ ,  is  significantly  less  than  the  power  saved  from  the 
drag  reduction,  A CflCfTwUc**3.2pu3ry  where  pw ,  p,  Cy, 
tw,  and  Uc  are  the  wall  pressure,  density,  skin-friction  co¬ 
efficient,  averaged  wall-shear  stress,  and  the  centerline  ve¬ 
locity,  respectively. 

Although  the  present  two-dimensional  controller  aug¬ 
mented  by  an  ad  hoc  three-dimensional  controller  has  shown 
a  promising  result,  it  is  apparent  that  we  need  to  develop  a 
three-dimensional  controller  using  the  same  formulation  pre¬ 
sented  in  this  paper.  Extensions  of  LQG(W2)/LTR  design  by 
using  three-dimensional  channel  flow  models  are  in 
progress. 
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